diff --git a/AculCalib/AculCalPars.cpp b/AculCalib/AculCalPars.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fcce0a249252c59591628e37c02fcd128e72a09f --- /dev/null +++ b/AculCalib/AculCalPars.cpp @@ -0,0 +1,40 @@ +/* + * AculCalParameters.cpp + * + * Created on: Oct 20, 2016 + * Author: vratik + */ + +#include "AculCalPars.h" + +AculCalPars::AculCalPars() { + // TODO Auto-generated constructor stub + Reset(); + +} + +AculCalPars::AculCalPars(const char* parFile) { + // TODO Auto-generated constructor stub + + SetParFile(parFile); + Init(); + +} +AculCalPars::~AculCalPars() { + // TODO Auto-generated destructor stub +} + +void AculCalPars::Init() { + SetPars(); +} + +void AculCalPars::SetParFile(const char* parfile) { + + fParFileName = parfile; + return; +} + +void AculCalPars::Reset() { + fParFileName = ""; + fEtab.clear(); +} diff --git a/AculCalib/AculCalPars.h b/AculCalib/AculCalPars.h new file mode 100644 index 0000000000000000000000000000000000000000..292343425722cbe83d06149fcc2821a568605e23 --- /dev/null +++ b/AculCalib/AculCalPars.h @@ -0,0 +1,74 @@ +/* + * AculCalParameters.h + * + * Created on: Oct 20, 2016 + * Author: vratik + */ + +#ifndef ACULCALIB_ACULCALPARS_H_ +#define ACULCALIB_ACULCALPARS_H_ + +#include +#include +#include + +#include "TString.h" +#include "TArrayD.h" +#include "TArrayI.h" +#include "TCutG.h" + +using std::cout; +using std::cerr; +using std::endl; +using std::vector; + +class AculCalPars { + +protected: +//general + TString fParFileName; + + vector fEtab; //energy which will be used for calibration + + //todo energies used for calibration should be in this class + + + +public: + AculCalPars(); + AculCalPars(const char* parFile); + virtual ~AculCalPars(); + ClassDef(AculCalPars,1); + + virtual void Init(); + + //getters + const char* GetParFileName() {return fParFileName.Data();} +// Int_t GetNoCrystals() {return fNoCrystals;} + virtual const char* GetDetName() {return 0;} + virtual const char* GetParticleName() {return 0;} + virtual Int_t GetNoRawFiles() {return 0;}; + virtual const char* GetFileName(Int_t i) {return 0;}; + virtual Int_t GetNoEPoints() {return 0;} + virtual Double_t GetCalEnergy(Int_t i) {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;}; + + void SetParFileName(const char* parFile) {fParFileName = parFile;} + + void SetParFile(const char* parfile); + virtual void PrintParameters(const char* option = "") {}; + + virtual void Reset(); + +protected: + virtual void SetPars() {}; +}; + +#endif /* ACULCALIB_ACULCALPARS_H_ */ diff --git a/AculCalib/AculCalParsScint.cpp b/AculCalib/AculCalParsScint.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fada257520dfc69b8146db7ef467b8d41567e0e9 --- /dev/null +++ b/AculCalib/AculCalParsScint.cpp @@ -0,0 +1,360 @@ +/* + * AculCalParsScint.cpp + * + * Created on: Oct 21, 2016 + * Author: vratik + */ + +#include "AculCalParsScint.h" +#include "TFile.h" + +AculCalParsScint::AculCalParsScint() { + // TODO Auto-generated constructor stub + Reset(); + +} + +AculCalParsScint::AculCalParsScint(const char* parFile) { + // TODO Auto-generated constructor stub + SetParFile(parFile); + Init(); + +} + +AculCalParsScint::~AculCalParsScint() { + // TODO Auto-generated destructor stub +} + +void AculCalParsScint::Init() { + SetPars(); +// LoadCuts(); +} + +void AculCalParsScint::SetPars() { + + std::ifstream infile(fParFileName.Data()); + if ( !infile.is_open() ) { + printf("AculCalParsScint::ReadParFile: File %s was not open and parameters were not set.\n", fParFileName.Data()); + Reset(); + return; + } + + fEnergyPoints = 0; + + TString line; + TString word; + + Int_t i, min, max; + Int_t lineLength = 400; + Char_t det[lineLength]; + Char_t part[lineLength]; + Char_t fname[lineLength]; + Char_t cname[lineLength]; + + double en; //energy + + while (!infile.eof()) { + line.ReadLine(infile); + + if ( line.BeginsWith("#") || line.BeginsWith("//") ) continue; + + if ( line.BeginsWith("energies") ) { + sscanf(line.Data(), "%*s %d", &i); + continue; + } + + if ( line.BeginsWith("crystals") ) { + sscanf(line.Data(), "%*s %d", &i); + fNoCrystals = i; + continue; + } + + if ( line.BeginsWith("files") ) { + sscanf(line.Data(), "%*s %d", &i); + fNoFiles = i; + fFilePars.resize(fNoFiles); + for (Int_t j = 0; j < fNoFiles; j++) { + line.ReadLine(infile); + if ( line.BeginsWith("#") || line.BeginsWith("//") ) { + j--; + continue; + } + sscanf(line, "%s", fname); + word = fname; + fFilePars[j].SetRawFileName(fname); + fFilePars[j].SetNoCrystals(fNoCrystals); + fFilePars[j].Init(); + } + continue; + } + + /*if ( line.BeginsWith("cutFile") ) { + sscanf(line.Data(), "%*s %s %d", fname, &fNoCuts); + fCutsFileName = fname; + for (Int_t j = 0; j < fNoCuts; j++) { + line.ReadLine(infile); + if ( line.BeginsWith("#") || line.BeginsWith("//") ) { + j--; + continue; + } + sscanf(line, "%s", cname); + fCutName.push_back(cname); + } + continue; + }*/ + + if ( line.BeginsWith("detector") ) { + sscanf(line.Data(), "%*s %s %s", det, part); + fDetName = det; + fPartName = part; + continue; + } + + 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); + + }//while + + infile.close(); + return; +} + +void AculCalParsScint::PrintParameters(const char* option) { + + TString opt = option; + + cout << "Parameters read from file \"" << fParFileName.Data() << "\"." << endl; + cout << "\tCalibration of detector \"" << fDetName << "\" with " << fNoCrystals << " crystals." << endl; + cout << "\tParticle: " << fPartName << endl; + cout << "\tInput files with raw data (" << fNoFiles << "):" << endl; + for (Int_t i = 0; i < (Int_t)fFilePars.size(); i++) { + cout << "\t " << GetFileName(i) << endl; + } + + cout << "\tInput energies (" << fNoFiles << "):" << endl; + 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; +// } +// } + + + + + 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; + } + } + + return; +} + +void AculCalParsScint::Reset() { + + fNoCrystals = 0; + + fParFileName = ""; + + fDetName = ""; + fPartName = ""; + + fNoFiles = 0; + fFilePars.clear(); + fEnergyPoints = 0; + +// fCutsFileName = ""; +// fNoCuts = 0; //number of cuts + +// fCutName.clear(); + + return; +} + +const char* AculCalParsScint::GetFileName(Int_t i) { + + if ( i >= (Int_t)fFilePars.size() ) { + cerr << "\"AculCalParsScint::GetFileName\" index i cannot be higher than " << fFilePars.size() - 1 << endl; + return 0; + } + return fFilePars[i].GetRawFileName(); +} + +/*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) { + 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, Int_t treeID) { + + const TString cName = cutName; + + 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; + } + } + + + cerr << "\"AculCalParsScint::GetCut\" cut \"" << cutName << "\" was not found." << endl; + return 0; +} + +Double_t AculCalParsScint::GetCalEnergy(Int_t i) { + + if ( i >= (Int_t)fFilePars.size() ) { + cerr << "\"AculCalParsScint::GetCalEnergy\" index i cannot be higher than " << fFilePars.size() - 1 << endl; + return 0; + } + return fFilePars[i].GetEnergy(); +} + +Int_t AculCalParsScint::GetMinChannel(Int_t energy, Int_t crystal) { + + if ( energy >= (Int_t)fFilePars.size() ) { + cerr << "\"AculCalParsScint::GetMinChannel\" index \"energy\" cannot be higher than " + << fFilePars.size() - 1 << endl; + return 0; + } + + if ( crystal >= fFilePars[energy].GetNoCrystals() ) { + cerr << "\"AculCalParsScint::GetMinChannel\" index \"crystal\" cannot be higher than " + << fFilePars[energy].GetNoCrystals() - 1 << endl; + return 0; + } + + return fFilePars[energy].GetPeakMin(crystal); +} + +Int_t AculCalParsScint::GetMaxChannel(Int_t energy, Int_t crystal) { + + if ( energy >= (Int_t)fFilePars.size() ) { + cerr << "\"AculCalParsScint::GetMinChannel\" index \"energy\" cannot be higher than " + << fFilePars.size() - 1 << endl; + return 0; + } + + if ( crystal >= fFilePars[energy].GetNoCrystals() ) { + cerr << "\"AculCalParsScint::GetMinChannel\" index \"crystal\" cannot be higher than " + << fFilePars[energy].GetNoCrystals() - 1 << endl; + return 0; + } + + return fFilePars[energy].GetPeakMax(crystal); +} + +/*void AculCalParsScint::LoadCuts() { + + if (fCutsFileName.Length() == 0) { + printf("AculCalParsScint::LoadCuts: Name of file (*.root) with cuts was not provided.\n"); + printf("AculCalParsScint::LoadCuts: No cuts has been loaded.\n"); + return; + } + + TFile cutFile(fCutsFileName.Data(), "READ"); + if(!cutFile.IsOpen()) { + cerr << "\"AculCalParsScint::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 << "\"AculCalParsScint::LoadCuts\" Cut \"" << fCutName[i] + << "\" was not found in file " << fCutsFileName << "." << endl; + continue; + } + fCuts.push_back(*currentCut); + } + +}*/ diff --git a/AculCalib/AculCalParsScint.h b/AculCalib/AculCalParsScint.h new file mode 100644 index 0000000000000000000000000000000000000000..90b7b5ea06b16a96b67450ea75ad2a6af4964bc8 --- /dev/null +++ b/AculCalib/AculCalParsScint.h @@ -0,0 +1,74 @@ +/* + * AculCalParsScint.h + * + * Created on: Oct 21, 2016 + * Author: vratik + */ + +#ifndef ACULCALIB_ACULCALPARSSCINT_H_ +#define ACULCALIB_ACULCALPARSSCINT_H_ + +#include "AculCalPars.h" +#include "AculCalParsScintFile.h" +//#include "TCutG.h" + +using std::cerr; + +class AculCalParsScint: public AculCalPars { + +private: + Int_t fNoCrystals; + + TString fDetName; + TString fPartName; + + Int_t fNoFiles; + 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; + + Int_t fEnergyPoints; + +// TString fCutsFileName; +// Int_t fNoCuts; //number of cuts + +public: + AculCalParsScint(); + AculCalParsScint(const char* parFile); + virtual ~AculCalParsScint(); + ClassDef(AculCalParsScint, 1); + + virtual void Init(); + + void PrintParameters(const char* option = ""); + // If option contains "all", all parameters will be printed out. + // By default important parameters only will be printed out. + + void Reset(); + + //getters + Int_t GetNoCrystals() {return fNoCrystals;} + const char* GetDetName() {return fDetName.Data();} + const char* GetParticleName() {return fPartName.Data();} + Int_t GetNoRawFiles() {return fNoFiles;} + const char* GetFileName(Int_t i); + 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, 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(); + +protected: + void SetPars(); +}; + +#endif /* ACULCALIB_ACULCALPARSSCINT_H_ */ diff --git a/AculCalib/AculCalParsScintFile.cpp b/AculCalib/AculCalParsScintFile.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cce8173bf3d64ac3f8ada7136ee1d2ae557bb734 --- /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 0000000000000000000000000000000000000000..db5ad8de4e53efcbd04907766e8e2ab780b1319b --- /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 0000000000000000000000000000000000000000..b580524fcdf1c85489ff2ddf0da059a701146c33 --- /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]; + + + std::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 0000000000000000000000000000000000000000..c4a0e6a071495e1e9d2718bbf614041225341b5d --- /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/AculCalib.cpp b/AculCalib/AculCalib.cpp new file mode 100644 index 0000000000000000000000000000000000000000..16b914744bd7eb83fbbd24f448e7ba38e81bcc5d --- /dev/null +++ b/AculCalib/AculCalib.cpp @@ -0,0 +1,96 @@ +#include "AculCalib.h" + +//---------------------------------------------------------------------------- +AculCalib::AculCalib() { +//documentation example + printf("AculCalib::Default constructor called.\n"); + + Init(); + +} + +//---------------------------------------------------------------------------- +/*AculCalib::AculCalib(const char* parfile) { + printf("AculCalib::Constructor called.\n"); + + SetParFile(parfile); + Init(); + +}*/ + +//---------------------------------------------------------------------------- +AculCalib::~AculCalib() { + + printf("AculCalib::Destructor called.\n"); + + delete fPars; + +} + +//---------------------------------------------------------------------------- +void AculCalib::Init() { + + fPars = 0; +} + +//---------------------------------------------------------------------------- +void AculCalib::SetParFile(const char* parfile) { + + fParFileName = parfile; + return; +} + +//---------------------------------------------------------------------------- +void AculCalib::PrintParameters(const char* option) { + if (!fPars) { + cerr << "\"AculCalib::PrintParameters\" parameters were not initialized." << endl; + return; + } + + fPars->PrintParameters(option); +} + +//---------------------------------------------------------------------------- +void AculCalib::PrintCalibParameters() { + cout << "i\ta\t\b" << endl; + for (Int_t i = 0; i<=(Int_t)fA.size(); i++) { + cout << i << "\t" << fA[i] << "\t" << fB[i] << endl; + } +} + +//---------------------------------------------------------------------------- +Double_t AculCalib::GetA(Int_t i) { + if (i >= (Int_t)fA.size()) //if i >= number of array element + { + return 0.; + } + return fA[i]; +} + +//---------------------------------------------------------------------------- +Double_t AculCalib::GetB(Int_t i) { + if (i >= (Int_t)fB.size()) //if i >= number of array element + { + return 0.; + } + return fB[i]; +} + +//---------------------------------------------------------------------------- +void AculCalib::CanDivider(TCanvas *c, Int_t noPads, Int_t columns, Int_t rows) { + + c->Clear(); + + if (noPads == 16) { + c->Divide(4, 4); + return; + } + + if (rows != 0 && columns != 0) { + c->Divide(columns, rows); + return; + } +} + +//---------------------------------------------------------------------------- +ClassImp(AculCalib); diff --git a/AculCalib/AculCalib.h b/AculCalib/AculCalib.h new file mode 100644 index 0000000000000000000000000000000000000000..782b5106b6d074fde1e06eb96c5cec7e43edf87e --- /dev/null +++ b/AculCalib/AculCalib.h @@ -0,0 +1,65 @@ +//////////////////////////////////////////////////////////////// +// // +// TMyClass // +// // +// This is the description block. // +// // +//////////////////////////////////////////////////////////////// + +#ifndef ACULCALIB_ACULCALIB_H_ +#define ACULCALIB_ACULCALIB_H_ + +//#include "TObject.h" +//#include "TROOT.h" +#include +#include + +#include "TArrayD.h" +#include "TString.h" +#include "TCanvas.h" + +#include "./AculCalPars.h" + +using std::cout; +using std::endl; +using std::cerr; + +class AculCalib { + +protected: + vector fA; //calibration parameter A + vector fB; //calibration parameter A + + TString fParFileName; //name of file with parameters + + AculCalPars *fPars; //parameters used for calibration + +public: + //Default constructor. + AculCalib(); + + //Destructor. + virtual ~AculCalib(); + ClassDef(AculCalib,1); + + virtual void Init(); + + void SetParFile(const char* parfile); + const char* GetParFileName() {return fParFileName.Data();} + + //Print parameters used for calibration. + void PrintParameters(const char* option = ""); + //Print calibration parameters. + void PrintCalibParameters(); + + //Get calibration parameter A. + Double_t GetA(Int_t i); + //Get calibration parameter B. + Double_t GetB(Int_t i); + +protected: + void CanDivider(TCanvas *c, Int_t noPads, Int_t columns = 0, Int_t rows = 0); + +}; + +#endif diff --git a/AculCalib/AculCalib.mk b/AculCalib/AculCalib.mk new file mode 100755 index 0000000000000000000000000000000000000000..4aab82ebd365326cd89349e38b8389fad461519a --- /dev/null +++ b/AculCalib/AculCalib.mk @@ -0,0 +1,47 @@ +################################################################################ +# AculData input with some variables +################################################################################ + +# ACULCALIBLIBS := -lCore -lCint -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf -lSpectrum #-lTELoss +ACULCALIBLIBS := -lCore -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf -lSpectrum #-lTELoss + +# Add inputs and outputs from these tool invocations to the build variables +ACULCALIB_HEADERS += \ +$(ACULCALIB)/AculCalib.h \ +$(ACULCALIB)/AculCalibScint.h \ +$(ACULCALIB)/AculCalibSi.h \ +$(ACULCALIB)/AculCalPars.h \ +$(ACULCALIB)/AculCalParsScint.h \ +$(ACULCALIB)/AculCalParsScintFile.h \ +$(ACULCALIB)/AculCalParsSi.h \ +$(ACULCALIB)/linkdef.h + +ACULCALIBCPP_SRCS += \ +$(ACULCALIB)/AculCalib.cpp \ +$(ACULCALIB)/AculCalibScint.cpp \ +$(ACULCALIB)/AculCalibSi.cpp \ +$(ACULCALIB)/AculCalPars.cpp \ +$(ACULCALIB)/AculCalParsScint.cpp \ +$(ACULCALIB)/AculCalParsScintFile.cpp \ +$(ACULCALIB)/AculCalParsSi.cpp \ +$(ACULCALIB)/AculCalibDict.cpp + +ACULCALIBOBJS += \ +$(ACULCALIB)/AculCalib.o \ +$(ACULCALIB)/AculCalibScint.o \ +$(ACULCALIB)/AculCalibSi.o \ +$(ACULCALIB)/AculCalPars.o \ +$(ACULCALIB)/AculCalParsScint.o \ +$(ACULCALIB)/AculCalParsScintFile.o \ +$(ACULCALIB)/AculCalParsSi.o \ +$(ACULCALIB)/AculCalibDict.o + +ACULCALIBCPP_DEPS += \ +$(ACULCALIB)/AculCalib.d \ +$(ACULCALIB)/AculCalibScint.d \ +$(ACULCALIB)/AculCalibSi.d \ +$(ACULCALIB)/AculCalPars.d \ +$(ACULCALIB)/AculCalParsScint.d \ +$(ACULCALIB)/AculCalParsScintFile.d \ +$(ACULCALIB)/AculCalParsSi.d \ +$(ACULCALIB)/AculCalibDict.d \ No newline at end of file diff --git a/AculCalib/AculCalibScint.cpp b/AculCalib/AculCalibScint.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0bed6fa47a490fbf2a9547a6b8a9239b4df00f96 --- /dev/null +++ b/AculCalib/AculCalibScint.cpp @@ -0,0 +1,426 @@ +#include "AculCalibScint.h" +#include "AculCalParsScint.h" + +ClassImp(AculCalibScint); + +AculCalibScint::AculCalibScint() { + + printf("AculCalibScint::Default constructor called.\n"); + fInFiles = 0; + fTrees = 0; + fCutFile = 0; + + +} + +AculCalibScint::AculCalibScint(const char* parfile) { + printf("AculCalibScint::Constructor called.\n"); + + SetParFile(parfile); + Init(); +// cout << "nofiles: " << nofiles << endl; +// OpenTrees(); +// LoadCuts(); + +// fr.Print(); +// fr.At(0); + + +} + +AculCalibScint::~AculCalibScint() { + + printf("AculCalibScint::Destructor called.\n"); + +} + +void AculCalibScint::Init() { + fPars = new AculCalParsScint(fParFileName.Data()); + OpenFiles(); + LoadTrees(); + + fHistFull.resize(fPars->GetNoRawFiles()); + fHistCut.resize(fPars->GetNoRawFiles()); + + fMeanPeakPos.resize(fPars->GetNoRawFiles()); + fMeanPeakRMS.resize(fPars->GetNoRawFiles()); +// SetPars(); +} + +void AculCalibScint::PrintTrees() { + + if (!fTrees) { + cerr << "\"AculCalibScint::PrintTrees\" Probably no tree was open." << endl; + return; + } + +// TTree *curTree = 0; + cout << " AculCalibScint::PrintTrees:" << endl; + for (Int_t i = 0; i < fPars->GetNoRawFiles(); i++) { +// curTree = fTrees[i]; + if (fTrees[i]) { + printf("\tTree No. %d; File: %s; Name: %s\n", i, fTrees[i]->GetDirectory()->GetName(), fTrees[i]->GetName()); + } else { + printf("\tTree No. %d was not loaded. Maximal number of trees is %d\n", i, fPars->GetNoRawFiles()); + } + } + + return; +} + +void AculCalibScint::OpenFiles() { + + if (!fPars) { + cerr << "\"AculCalibScint::OpenFiles\" parameters were not initialized." << endl; + return; + } + + fPars->GetNoRawFiles(); + fInFiles = new TFile* [fPars->GetNoRawFiles()]; + for (Int_t i = 0; i < fPars->GetNoRawFiles(); i++) { + fInFiles[i] = new TFile(fPars->GetFileName(i), "READ"); +// cout << "\"AculCalibScint::OpenFiles\" File \"" << fInFiles[i]->GetName() << "\" was opened." << endl; + } + +} + +void AculCalibScint::LoadTrees() { + + if (!fInFiles) { + cerr << "\"AculCalibScint::LoadTrees\" Input files were not open." << endl; + return; + } + + fTrees = new TTree* [fPars->GetNoRawFiles()]; + for (Int_t i = 0; i < fPars->GetNoRawFiles(); i++) { + fTrees[i] = (TTree*)fInFiles[i]->Get("AnalysisxTree"); +// cout << "\"AculCalibScint::LoadTrees\" Tree \"" << fTrees[i]->GetName() +// << "\" from file \"" << fInFiles[i]->GetName() << "\" was loaded." << endl; + } + +} + +TTree* AculCalibScint::GetTree(Int_t treeID) { + if (treeID >= fPars->GetNoRawFiles()) { + cerr << "\"AculCalibScint::DrawVariable\" Tree number " << treeID + << " cannot exist. Maximal treeID is " << fPars->GetNoRawFiles() - 1 << "." << endl; + return 0; + } + + return fTrees[treeID]; +} + +void AculCalibScint::DrawVariable(const char* variable, Int_t treeID, TCanvas *canvas, Int_t lowRange, Int_t upRange) { + + + if (!canvas) { + cerr << "\"AculCalibScint::DrawVariable\" Canvas pointer was NULL" << endl; + return; + } + + //todo change 16 to parameter + const Int_t nPads = 16; + CanDivider(canvas, nPads); + + TString canvasTitle; + TString var; + TString con; + + TTree *curTree = 0; + curTree = GetTree(treeID); + if (!curTree) { + printf("AculCalibScint::DrawVariable: Tree No. %d was not found.\n", treeID); + return; + } + + canvasTitle.Form("variable: %s; tree: %d", variable, treeID); + canvas->SetTitle(canvasTitle.Data()); + + for (Int_t i = 0; i < nPads; i++) { + var.Form("%s[%d]", variable, i); + con.Form("%s[%d]>%d && %s[%d]<%d", variable, i, lowRange, variable, i, upRange); + canvas->cd(i+1); + curTree->Draw(var.Data(), con.Data()); + canvas->Update(); + } + +} + +void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable) { + + canvas->SetTitle("Beam"); + + const Int_t noFiles = file; + + //todo delet parameter "files" + CanDivider(canvas, 0, noFiles, 3); + + TTree *curTree = 0; + TCutG *curCutG = 0; + TString var; + TString con; + +// 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); + curTree = GetTree(i); + if (!curTree) { + printf("AculCalibCsI::DrawBeam: Tree No. %d was not found.\n", i); + continue; + } + curTree->Draw("QDC[0]:TDC[0]", "TDC[0]<1000 && QDC[0]<2000", "cont"); +// 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); + 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); + 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); + var.Form("%s[5]:TDC[0]", variable); + con.Form("%s[5]>200", variable); + curTree->Draw(var.Data(), con.Data(), "cont"); + for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) { +// 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(); + } +} + +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; + CanDivider(canvas, noCrystals); + + TString var; + 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 + 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); + if (!curTree) { + printf("AculCalibCsI::GetPeakMean: Tree No. %d was not found.\n", treeID); + return; + } + + TString sVariable = variable; + sVariable.ToLower(); + Int_t channel = 0; + + for (Int_t i = 0; i < noCrystals; i++) { + + if (sVariable.Contains("anc")) { + channel = i+1; + } else { + channel = i; + } + + var.Form("%s[%d]>>hfull[%d][%d]", variable, channel, treeID, i); + 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) ); + curTree->SetLineColor(1); + curTree->Draw(var.Data(), con.Data()); + + var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, treeID, i); + 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); + curTree->Draw(var.Data(), con.Data(), "same"); + + gPad->Update(); + fMeanPeakPos[treeID].push_back(fHistCut[treeID][i]->GetMean()); +// cout << fMeanPeakPos[treeID][i] << endl; + fMeanPeakRMS[treeID].push_back(fHistCut[treeID][i]->GetRMS()); + } + + +// cout << fMeanPeakPos.size() << "\t" << fMeanPeakPos[treeID].size() << endl; +// cout << fMeanPeakPos.size() << "\t" << fMeanPeakPos[0].size() << endl; +// cout << fMeanPeakPos.size() << "\t" << fMeanPeakPos[1].size() << endl; +// cout << fMeanPeakPos.size() << "\t" << fMeanPeakPos[2].size() << endl; +// cout << fMeanPeakPos.size() << "\t" << fMeanPeakPos[3].size() << endl; + + canvas->Update(); + +} + +void AculCalibScint::Calibrate(TCanvas *canvas, Bool_t savefile, const char* filename, const char* option) { + + //todo change this parameter as nonconstant + //e.g. the number of found fPeakMeans + const Int_t noChannels = 16; + + //todo + if (!canvas) return; + CanDivider(canvas, noChannels); + + for (Int_t i = 0; i < fPars->GetNoRawFiles(); i++) { + cout << i << "\t" << fMeanPeakPos[i].size() << endl; + if (fMeanPeakPos[i].size() < noChannels || fMeanPeakRMS[i].size() < noChannels) { +// cout << i << "\t" << fMeanPeakPos[i].size() << endl; + cout << "\"AculCalibScint::Calibrate\" something bad here" << endl; + return; + } + } + +// cout << alphas2.GetSize()+1 << endl; +// cout << energyPoints+1 << endl; + cout << fPars->GetNoEPoints() << endl; + const Int_t energyPoints = fPars->GetNoEPoints(); + + TString gName; + TString gTitle; + +// if (savefile) fGraphs->Open(filename, "RECREATE"); + + TF1 *fnc; + + for (Int_t i = 0; i < noChannels; i++) { + canvas->cd(i+1); + gCal[i] = new TGraphErrors(energyPoints+1); +// FillGraph(gCal[i], energies.GetSize()+1, energies.GetArray(), i); + FillGraph(gCal[i], energyPoints+1, i); +// if (savefile) gCal[i]->Write(); +// gCal[i]->Draw("Al*"); + gCal[i]->Draw("A*"); + gName.Form("g%s%s%d\n", fPars->GetDetName(), fPars->GetParticleName(), i); + gTitle.Form("%s %s\n", fPars->GetDetName(), fPars->GetParticleName()); +// gCal[i]->SetTitle(gTitle.Data()); + gCal[i]->SetName(gName.Data()); + gCal[i]->Fit("pol1"); + fnc = gCal[i]->GetFunction("pol1"); + fnc->SetLineColor(kRed); + fA.push_back(fnc->GetParameter(1)); + fB.push_back(fnc->GetParameter(0)); + canvas->Update(); + } + if (savefile) SaveClbGraphs(filename, option); + +} + +void AculCalibScint::FillGraph(TGraphErrors *g, Int_t npoints, /*Double_t *energies,*/ Int_t graphNumber, const char* option) { + + TString opt = option; + + //todo check if fMeanPeakPos is full + + //all available energy points and (0,0) + g->SetPoint(0, 0., 0.); + 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]); +// cout << "Graph number: " << graphNumber << endl; +// cout << graphNumber << "\t" << fMeanPeakPos[i][graphNumber] << endl; + g->SetPoint(i+1, fMeanPeakPos[i][graphNumber], fPars->GetCalEnergy(i)); + g->SetPointError(i+1, fMeanPeakRMS[i][graphNumber], 0); + } + +// for (Int_t j = 1; j < 4; j++) { +// g->SetPoint(j, 0., 0.); +// } + + +} + +void AculCalibScint::SaveClbGraphs(const char* filename, const char* option) { + +// cout << "asdasd" << endl; +// cout << fGraphs << endl; + TFile fGraphs(filename, option); + +// if (fGraphs) fGraphs->Close(); + +// cout << "asdasd" << endl; +// fGraphs = new TFile(filename, option); + cout << fGraphs.IsOpen() << endl; + cout << fGraphs.GetName() << endl; + fGraphs.Print(); +// if (!fGraphs->IsOpen()) { +// printf("AculCalibCsI::SaveClbGraphs: file %s was not open.\n", filename); +// return; +// } + + for (Int_t i = 0; i<16; i++) { + fGraphs.cd(); + gCal[i]->Write(); + } + fGraphs.Close(); + return; +} + +void AculCalibScint::WriteClbParameters(const char* filename) { + + std::ofstream outfile(filename); + if ( !outfile.is_open() ) { + printf("AculCalibCsI::WriteClbParameters: File %s was not open.\n", filename); + return; + } + + outfile << "#detector:\t" << fPars->GetDetName() << ",\tparticle:\t" << fPars->GetParticleName() << endl; + outfile << "#channel\tfA\tfB" << endl; + for (Int_t i = 0; i < (Int_t)fA.size(); i++) { + outfile << i << "\t" << fA[i] << "\t" << fB[i] << endl; + } + + outfile.close(); +} diff --git a/AculCalib/AculCalibScint.h b/AculCalib/AculCalibScint.h new file mode 100644 index 0000000000000000000000000000000000000000..4e3682221de9a166b5430d50ddeb4ca0a7630320 --- /dev/null +++ b/AculCalib/AculCalibScint.h @@ -0,0 +1,134 @@ +////////////////////////////////////////////////////////////////////////////// +// // +// This class is intended for ... // +// // +// It work like this... User may use such functions ... in such steps. // +// +// Keywords in parameter file and possible comments +// +// Something about order of keywords +// +// Example of parameter file (with explanations) +// +// energies 4 #number of energy points +// crystals 16 #number of crystals in detector +// files 4 +// #fileID fileName +// #fileID not implemented yet +// data/csi_13_Ealpha16.root +// data/csi_13_Ealpha21.root +// data/csi_13_Ealpha26.root +// data/csi_13_Ealpha30.root +// #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 +//\ ... +// 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 +//\ \... +// 15 1150 1350 +// #CsI2, energy ~ 26 MeV/A, Alpha +// energy 105.5 +// ... +// +//asdasd +////////////////////////////////////////////////////////////////////////////// + + +#pragma once + +//#include "TObject.h" +//#include "TROOT.h" +#include +#include +//#include + +#include "TFile.h" +#include "TTree.h" +//#include "TCanvas.h" +#include "TH1I.h" +#include "TGraphErrors.h" +//#include "TArrayD.h" +#include "TF1.h" + +#include "./AculCalib.h" + +using std::cout; +using std::endl; +//using std::vector; + +class AculCalibScint : public AculCalib { + +private: +// + //todo delete this strange double array somewhere + TFile **fInFiles; //! doc of class member + TTree **fTrees; //! + + TFile *fCutFile; + + vector< vector > fHistFull; + vector< vector > fHistCut; +// + vector< vector > fMeanPeakPos; //doc of class member + vector< vector > fMeanPeakRMS; +// + //todo make vector of graphs + TGraphErrors *gCal[16]; +// TFile *fGraphs; + +public: + AculCalibScint(); + //default constructor + + AculCalibScint(const char* parfile); + virtual ~AculCalibScint(); + ClassDef (AculCalibScint,1); //decription right after ClassDef + + virtual void Init(); + + void PrintTrees(); + + //Description of the function. I think it does not work. + void DrawVariable(const char* variable, Int_t treeID, TCanvas *canvas, Int_t lowRange = 0, Int_t upRange = 4096); + + //Description of the function. I think it does not work. + void DrawBeam(TCanvas *canvas, Int_t file, const char* variable); + + //One of the crucial function. It founds the mean of the peak related to our particle. + void GetPeakMean(const char* variable, Int_t treeID, 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); + + TTree* GetTree(Int_t treeID); + +//private functions: +// void LoadCuts(); + +private: + void OpenFiles(); + void LoadTrees(); + + void SaveClbGraphs(const char* filename, const char* option = "READ"); + void FillGraph(TGraphErrors *g, Int_t npoints, /*Double_t *energies,*/ Int_t graphNumber, const char* option = ""); + + +// void SetPars(); +}; diff --git a/AculCalib/AculCalibSi.cpp b/AculCalib/AculCalibSi.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0ffbc6f6a0681d28007695ba76d76a2d70159d57 --- /dev/null +++ b/AculCalib/AculCalibSi.cpp @@ -0,0 +1,35 @@ +/* + * AculCalibSi.cpp + * + * Created on: Oct 26, 2016 + * Author: vratik + */ + +#include "AculCalibSi.h" +#include "AculCalParsSi.h" + +ClassImp(AculCalibSi); + +AculCalibSi::AculCalibSi() { + printf("AculCalibScint::Default constructor called.\n"); +// fInFiles = 0; +// fTrees = 0; +// fCutFile = 0; +} + +AculCalibSi::AculCalibSi(const char* parfile) { + printf("AculCalibScint::Constructor called.\n"); + + SetParFile(parfile); + Init(); +} + +AculCalibSi::~AculCalibSi() { + + printf("AculCalibScint::Destructor called.\n"); + +} + +void AculCalibSi::Init() { + fPars = new AculCalParsSi(fParFileName.Data()); +} diff --git a/AculCalib/AculCalibSi.h b/AculCalib/AculCalibSi.h new file mode 100644 index 0000000000000000000000000000000000000000..7db4e3e3c763bed4f566784459f1922cf2daf06c --- /dev/null +++ b/AculCalib/AculCalibSi.h @@ -0,0 +1,44 @@ +/* + * AculCalibSi.h + * + * Created on: Oct 26, 2016 + * Author: vratik + */ + +#ifndef ACULCALIB_ACULCALIBSI_H_ +#define ACULCALIB_ACULCALIBSI_H_ + +#include "AculCalib.h" + +class AculCalibSi: public AculCalib { +private: + +// Double_t something; + +public: + AculCalibSi(); + AculCalibSi(const char* parfile); + virtual ~AculCalibSi(); + ClassDef(AculCalibSi, 1); + + void Init(); + +// Bool_t CalculateCalibParameters(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); //calculate calibration parameters for given block in given file + //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 + + +}; + +#endif /* ACULCALIB_ACULCALIBSI_H_ */ diff --git a/AculCalib/linkdef.h b/AculCalib/linkdef.h new file mode 100755 index 0000000000000000000000000000000000000000..b444014dde9dbceef9bc0c0310cf8ed1bbdacd8e --- /dev/null +++ b/AculCalib/linkdef.h @@ -0,0 +1,17 @@ +#ifdef __CINT__ +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +//#pragma link C++ class AculRaw; +//#pragma link C++ class AculConvert; +#pragma link C++ class AculCalib; +#pragma link C++ class AculCalibScint; +#pragma link C++ class AculCalibSi; +#pragma link C++ class AculCalPars; +#pragma link C++ class AculCalParsScint; +#pragma link C++ class AculCalParsScintFile; +#pragma link C++ class AculCalParsSi; + +#endif + diff --git a/README b/README index 9e04925efe4495953f34e0208d6311ac8a4176bb..ee77b5340a77c18e1278bd18f996130e91d25b97 100755 --- a/README +++ b/README @@ -1,3 +1,6 @@ +Tested with +root_6.28_04 (precompiled binaries for Ubuntu 20.04) + I. Use with GNU make 1) set path in makefile according to your operating system diff --git a/Utilities/CalPars.cpp b/Utilities/CalPars.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5d30da42fe891633593146cca367325e8882422c --- /dev/null +++ b/Utilities/CalPars.cpp @@ -0,0 +1,93 @@ +#include "CalPars.h" + +ClassImp(CalPars); + + +CalPars::CalPars(const char* clbfile) { + + ReadParFile(clbfile); + +} + +CalPars::~CalPars() { + + printf("CalPars::Destructor called.\n"); + +} + +Double_t CalPars::GetA(Int_t i) { + + if (i >= fA.GetSize()) { + return 0.; + } + return fA[i]; +} + +Double_t CalPars::GetB(Int_t i) { + + if (i >= fB.GetSize()) { + return 0.; + } + return fB[i]; +} + +Double_t CalPars::GetC(Int_t i) { + + if (i >= fC.GetSize()) { + return 0.; + } + return fC[i]; +} + +void CalPars::ReadParFile(const char* clbfile) { + + clbFile = clbfile; + + std::ifstream infile(clbfile); + if ( !infile.is_open() ) { + printf("CalPars::ReadParFile: File %s was not open.\n", clbfile); + printf("CalPars::ReadParFile: Calibration parameters are empty!!!\n"); + return; + } + + TString line; + Int_t nopars; + Double_t a, b; + Int_t c; //threshold + + printf("CalPars::ReadParFile: File %s was open.\n", clbfile); + + line.ReadLine(infile); + sscanf(line.Data(), "%*d %d", &nopars); +// Info("CalPars::ReadParFile", "%d calibration parameters will be loaded", nopars); + printf("CalPars::ReadParFile: %d calibration parameters will be loaded\n", nopars); + + fA.Set(nopars); + fB.Set(nopars); + fC.Set(nopars); + + for (Int_t i = 0; i +#include + +#include "TArrayD.h" +#include "TArrayI.h" +#include "TString.h" +#include "TObject.h" + +using std::cout; +using std::endl; + +class CalPars /*: TObject*/ { + +private: + TArrayD fA; + TArrayD fB; + TArrayI fC; + + TString clbFile; + +public: + + CalPars() : fA(0), fB(0), fC(0) {}; + CalPars(const char* clbfile); + virtual ~CalPars(); + + Double_t GetA(Int_t i); + Double_t GetB(Int_t i); + Double_t GetC(Int_t i); + + void ReadParFile(const char* clbfile); + void PrintClbPars(); + + ClassDef(CalPars, 1); +}; diff --git a/Utilities/ConfigDictionary.cpp b/Utilities/ConfigDictionary.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7e3790be97eed63747cb27fd329887a292fd37fa --- /dev/null +++ b/Utilities/ConfigDictionary.cpp @@ -0,0 +1,229 @@ +#include "ConfigDictionary.h" + +ClassImp(ConfigDictionary); + +////////////////////////////////////////////////////////////////////////////// +// BEGIN_HTML +//

Config(uration)Dictionary class

+//
+// Author: Bartlomiej Hnatio 2012-08-06 +//

+// This is very useful class to convert strings containing pairs "key"="value" +// into fast dictionaries. This strings are often read from external files. +// From dictionary created in such way one can easily +// extract values in few supported formats (integer, double, boolean, string) +// It can also be used in opposite direction: when dictionary is created with +// given values and keys, one can create config string, which is most +// convenient to write in external files. +//


+// Most simple example of usage: +//

Suppose you have two variables fD and fI: +//
-------------------
+//    Double_t fD;
+//    Int_t fI;
+//
+//fD = 3.14;
+//fI = 2012;
+//-----------------------
+// To save its parameters into file you can create ConfigDictionary class +// instance and use SetDouble and SetInt functions to +// insert parameters values with arbitrarly defined keys (let them be +// really long and descriptive in this example): +//
---------------------
+//ConfigDictionary CD;
+//CD.SetDouble("Most important variable D",fD);
+//CD.SetInt("Equally important integer I",fI);
+//---------------------
+// Now configuration string saved in ConfigDictionary CD +// can be obtained with ToString method and should look like: +//
+//"Most important variable D"="3.14" "Equally important integer I"="2012"
+//
+//
It can be easily saved to a file using simplest fstream +// methods. And the advantage is that as a key one can use any string, +// so configuration file created in such way is very self-explanatory. +//
+//
+// Now lets suppose the opposite action - loading config from file. +// Imagine, that you have 1000 objects of A class, which config was saved +// to file - one object per line: +//
-----------------
+//"Most important variable D"="3.14" "Equally important integer I"="2012"
+//"Equally important integer I"="1011" "Most important variable D"="8.15"
+//"Most important variable D"="13.16" "Equally important integer I"="10"
+//(...)
+//----------------
+// Please notice that order in which pairs of keys and values are placed +// in file doesn't make any difference to the ConfigDictionary class. +// To recreate objects, you just read each line containing single +// config string and then create with it ConfigDictionary object +// using special constructor with string as argument, +// or using FromString method. After this, you can +// extract parameters (using GetDouble and GetInt methods) +// with using same keys as were used for saving, and ConfigDictionary +// will return their values: +//
----------------------
+//string line;//This line you read from file
+//ConfigDictionary CD(line);
+//fD = CD.GetDouble("Most important variable D");
+//fI = CD.GetInt("Equally important integer I");
+//--------------------
+//
+// And last but not least: what happens, if key requested by user doesn't +// exist in dictionary? (This can be caused by many reasons, mostly errors on +// user side, but not always). +// When you try to extract non-existent key using Get functions, +// an exception is risen. In normal program it should end execution and print +// some information about where program stopped working. +// But lets say that you don't want that, i.e. program may use default +// configs instead of those from files, or not all keys were that important. +// You can surround all uses of Get methods with try/catch clause. +// So when exception is risen, you will catch it, and decide, if you want +// the program to stop running, or anything else. Simple example is +// shown below: +//
-------------------
+//ConfigDictionary CD(some_string);
+//try{//try to read important variables:
+//  double d = CD.GetDouble("crucial var");
+//  bool b = CD.GetBool("most important b");
+//  int i = CD.GetInt("unique ID");
+//  //...and so on
+//}catch(std::string & e){//catch any kind of exception
+//  Error("Some crucial variable wasn't read, ending program!");
+//  return SOME_REALLY_BAD_ERROR_CODE;
+//}
+//try{//now less important, or optional:
+//  string name = CD.GetString("least important variable");
+//  double p = CD.GetDouble("optional parameter");
+//  //...and so on
+//}catch(std::string & f){
+//  Info("Some optional variables wasn't read!");
+//}
+//--------------------
+// +// ADDED 21.03.2012: +// +// With use of templates, now ConfigDictionary can handle complex objects +// reading and writing instantly as long as proper >> and << operators are +// provided for an object initialization. +// +// END_HTML +////////////////////////////////////////////////////////////////////////////// + +//_____________________________________________________________________________ +ConfigDictionary::ConfigDictionary(){ + //Empty map... +} + +//_____________________________________________________________________________ +ConfigDictionary::ConfigDictionary(std::string params){ + //Creates dictionary using FromString method + FromString(params); +} + +//_____________________________________________________________________________ +std::string ConfigDictionary::ToString(){ + //Builds string that can be easily saved to file. + //Same format that can be read from files or from QtGui + //This should work same way for all uses + + std::map::iterator it;//Iterator to map elements + std::stringstream ss;//stream helpful with adding strings + //iterate whole dictionary: + for (it=configMap.begin();it!=configMap.end();++it){ + //insert pairs "key"="value": + ss<<"\""<first<<"\""<<"="<<"\""<second<<"\""<<" "; + } + return ss.str(); +} + +//_____________________________________________________________________________ +void ConfigDictionary::FromString(std::string params){ + //params - TString containing list of key=value pairs + // + //Changes string formatted: + //"key1"="value1" "key2"="value with spaces 2" ... + //into map with keys and values + //Useful in lots of I/O functions, + //when we want to have a nice, readable format of data + + std::stringstream loading(params); + std::string k,v; + + while(!loading.fail()){ + getline(loading,k,'\"');//removes everything to first " + getline(loading,k,'\"');//All chars between first two "" are the key + getline(loading,v,'\"');//removes all until third " + getline(loading,v,'\"');//All between another pair of "" is the value + if (!loading.fail()) + configMap[k]=v; + } +} + +//_____________________________________________________________________________ +std::string ConfigDictionary::GetString(std::string key)throw(std::string){ + //Extracts string from given key + //(if it exist, otherwise raises exception) + //This works also for long strings with white spaces in it! + if (configMap.find(key) == configMap.end()){ + Error("ConfigDictionary::GetString", + "Couldn't find the key: %s!",key.c_str()); + throw(key); + } + return configMap[key]; +} + +//_____________________________________________________________________________ +int ConfigDictionary::GetInt(std::string key)throw(std::string){ + //Extracts integer from given key + //(if it exist, otherwise raises exception) + return Get(key); +} + +//_____________________________________________________________________________ +double ConfigDictionary::GetDouble(std::string key)throw(std::string){ + //Extracts double precision floating number from given key + //(if it exist, otherwise raises exception) + return Get(key); +} + +//_____________________________________________________________________________ +bool ConfigDictionary::GetBool(std::string key)throw(std::string){ + //Extracts boolean from given key + //(if it exist, otherwise raises exception) + //Differs from template to make it more readable. + if (configMap.find(key) == configMap.end()){ + Error("ConfigDictionary::GetBool", + "Couldn't find the key: %s!",key.c_str()); + throw(key); + } + //Convert string to bool: + if (configMap[key].compare("true") == 0) + return true; + else return false; +} + +//_____________________________________________________________________________ +void ConfigDictionary::SetString(std::string key,std::string value){ + //Sets value to key, no comments needed here... + configMap[key] = value; +} + +//_____________________________________________________________________________ +void ConfigDictionary::SetDouble(std::string key,double value){ + Set(key,value); +} + +//_____________________________________________________________________________ +void ConfigDictionary::SetInt(std::string key,int value){ + //Sets value to key, converts int to string first + Set(key,value); +} + +//_____________________________________________________________________________ +void ConfigDictionary::SetBool(std::string key,bool value){ + //Sets value to key, converts bool to string first + //Differs from template to make it more readable. + if (value == true) configMap[key] = "true"; + else configMap[key] = "false"; +} diff --git a/Utilities/ConfigDictionary.h b/Utilities/ConfigDictionary.h new file mode 100644 index 0000000000000000000000000000000000000000..a6fd7d457eb590b03ee6f3bb2d11f9d045c4c44c --- /dev/null +++ b/Utilities/ConfigDictionary.h @@ -0,0 +1,128 @@ +#ifndef CONFIGDICTIONARY_H +#define CONFIGDICTIONARY_H + + +#include "TObject.h" +#include "ReturnCodes.h" +#include "TStopwatch.h" +#include "TROOT.h" +#include +#include +#include +#include +#include "TLorentzVector.h" +#include +#include "TMath.h" +#include "TString.h" +#include "TTree.h" + + +#define ARRAY_DELIM '$' + +class ConfigDictionary{ +public: + typedef std::map::iterator CDIter; + ConfigDictionary(); + ConfigDictionary(std::string); + virtual ~ConfigDictionary(){};//empty virtual destructor + + ClassDef(ConfigDictionary,1); + std::string ToString(); + void FromString(std::string); + + //These throw errors if couldn't find key: + //Most of methods below are just using templates to + //make them usable in Cint. + std::string GetString(std::string)throw(std::string); + int GetInt(std::string)throw(std::string); + double GetDouble(std::string)throw(std::string); + bool GetBool(std::string)throw(std::string); + + //These will always set 'something' into map: + void SetString(std::string,std::string); + void SetDouble(std::string,double); + void SetInt(std::string,int); + void SetBool(std::string,bool); + + CDIter Begin(){return configMap.begin();}; + CDIter End(){return configMap.end();}; + + + //_________________________________________________________________________ + template + T Get(std::string key) throw (std::string) { + //This method will get object of any class + //that was - presumably - stored before in string format in CD. + //Object needs to have proper >> operator. + //If operation fails, returned item should have default value + //depending on its empty constructor. + //For strings containing spaces it won't work properly, as it + //would extract only first space followed item! + if (configMap.find(key) == configMap.end()) { + Error( "ConfigDictionary::Get", + "Couldn't find the key: %s!", + key.c_str()); + throw(key); + } + std::istringstream inStream(configMap[key]); + T item; + inStream >> item; + return item; + } + + + //_____________________________________________________________________________ + template + void Set(std::string key,T value){ + //Sets value to key, converts T to string first with << operator. + std::stringstream ss; + ss< + int GetArray(std::string key,std::vector & output)throw(std::string){ + //Preconditions: T must have >> operator + //Effect: searches for key in map and pushes values stored in it into output vector + //Returns: number of elements pushed into vector + int retNumber = 0; + if (configMap.find(key) == configMap.end()){ + Error( "ConfigDictionary::GetInt", + "Couldn't find the key: %s!",key.c_str()); + throw(key); + } + std::istringstream arrayString(configMap[key]); + std::string itemString; + //Dangerous - maybe? - exceptions for ss not handled... + while (arrayString){ + std::getline(arrayString,itemString,ARRAY_DELIM); + std::stringstream itemStream(itemString); + if (itemString.length() > 0){ + T item; + itemStream>>item; + output.push_back(item); + retNumber++; + } + } + return retNumber; + } + + //_________________________________________________________________________ + template + void SetArray(std::string key,std::vector & input){ + std::ostringstream arrayString; + for (unsigned int i=0;i configMap; +}; + +#endif diff --git a/Utilities/CsICalib.cpp b/Utilities/CsICalib.cpp new file mode 100644 index 0000000000000000000000000000000000000000..35c36d9e7ff3b3cac727dee39b3919bf644041ff --- /dev/null +++ b/Utilities/CsICalib.cpp @@ -0,0 +1,633 @@ +#include "CsICalib.h" + +ClassImp(CsICalib); + +CsICalib::CsICalib() : fr("TFile", 5), tr("TTree", 5), cutsCol("TCutG"), fA(16), fB(16) { + + printf("CsICalib::Default constructor called.\n"); + + +} + +CsICalib::CsICalib(const char* parfile) : fr("TFile", 5), tr("TTree", 5), cutsCol("TCutG", 10), fA(16), fB(16) { + printf("CsICalib::Constructor called.\n"); + + ReadParFile(parfile); +// PrintParameters(); + cout << "nofiles: " << nofiles << endl; + OpenTrees(); + LoadCuts(); + +// fr.Print(); +// fr.At(0); + + +} + +CsICalib::~CsICalib() { + + printf("CsICalib::Destructor called.\n"); + +} + +void CsICalib::OpenTrees() { + + TFile *file; + for (Int_t i = 0; i < nofiles; i++) { + fr[i] = new TFile(fileNames[i], "READ"); + file = (TFile*)fr.At(i); + cout << file->GetName() << endl; + tr[i] = (TTree*)file->Get("AnalysisxTree"); + } +} + +void CsICalib::LoadCuts() { + + if (cutsFileName == "") { + printf("CsICalib::LoadCuts: Name of file (*.root) with cuts was not provided.\n"); + printf("CsICalib::LoadCuts: No cuts has been loaded.\n"); + return; + } + fCuts = new TFile(cutsFileName.Data(), "READ"); + for (Int_t i = 0; i < nCuts; i++) { + cutsCol[i] = (TCutG*)fCuts->Get(cutNames[i]); + } + +} + +void CsICalib::ReadParFile(const char* parfile) { + + std::ifstream infile(parfile); + if ( !infile.is_open() ) { + printf("CsICalib::ReadParFile: File %s was not open.\n", parfile); + return; + } + + energyPoints = 0; + + TString line; + + Int_t i, min, max; + char det[100], part[100], fname[500], cname[100]; + + double en; //energy + + while (!infile.eof()) { + line.ReadLine(infile); + + if ( line.BeginsWith("#") ) continue; + + if ( line.BeginsWith("energies") ) { + sscanf(line.Data(), "%*s %d", &i); + fE.Set(i); + continue; + } + + if ( line.BeginsWith("files") ) { + sscanf(line.Data(), "%*s %d", &i); + nofiles = i; + printf("CsICalib::ReadParFile: %d files will be loaded:\n", nofiles); + for (Int_t j = 0; j < nofiles; j++) { + line.ReadLine(infile); + sscanf(line, "%s", fname); + fileNames[j] = fname; + cout << fileNames[j] << endl; + } + continue; + } + + if ( line.BeginsWith("cutFile") ) { + sscanf(line.Data(), "%*s %s %d", fname, &nCuts); + cutsFileName = fname; + for (Int_t j = 0; j < nCuts; j++) { + line.ReadLine(infile); + sscanf(line, "%s", cname); + cutNames[j] = cname; + cout << cutNames[j] << endl; + } + continue; + } + + if ( line.BeginsWith("detector") ) { + sscanf(line.Data(), "%*s %s %s", det, part); + detName = det; + partName = part; + printf("%s %s\n", detName.Data(), partName.Data()); + continue; + } + + if ( line.BeginsWith("energy") ) { + sscanf(line.Data(), "%*s %lf", &en); + fE[energyPoints] = en; + printf("%f\n", fE[energyPoints]); + energyPoints++; + continue; + } + + sscanf(line.Data(), "%d %d %d", &i, &min, &max); + peakMin[energyPoints-1][i] = min; + peakMax[energyPoints-1][i] = max; +// printf("%d %d %d\n", i, peakMin[energyPoints-1][i], peakMax[energyPoints-1][i]); + }//while + + printf("energyPoints: %d\n", energyPoints); + + infile.close(); + return; +} + +void CsICalib::PrintParameters(const char* option) { + + TString opt = option; + + printf("Energy points: %d\n", energyPoints); +// printf("Detector: %s, Particle: %s\n", ????); + for (Int_t i = 0; i < energyPoints; i++) { + printf("Peak ranges for energy %f:\n", fE[i]); + if ( opt.Contains("v") ) { + for (Int_t j = 0; j < 16; j++) { + printf("ch: %d\tmin: %d\tmax: %d\n", j, peakMin[i][j], peakMax[i][j]); + } + }//if + }//for + +} + + +void CsICalib::PrintPeakRanges() { + for (Int_t i = 0; i < 16; i++) { + printf("%d\t%d\t%d\n", i, peakMin[0][i], peakMax[0][i]); + } +} + +void CsICalib::PrintTrees() { + + TTree *curTree = 0; + for (Int_t i = 0; i < tr.GetEntries(); i++) { + curTree = (TTree*)tr.At(i); + if (curTree) { + printf("Tree No. %d; File: %s; Name: %s\n", i, curTree->GetDirectory()->GetName(), curTree->GetName()); + } else { + printf("Tree No. %d was not loaded. Maximal number of trees is %d\n", i, NOCALFILES); + } + } + + return; +} + +void CsICalib::PrintFiles() { + + printf("Number of loaded files: %d\n", fr.GetEntries()); + + + TFile *curFile = 0; + for (Int_t i = 0; i < fr.GetEntries(); i++) { + curFile = (TFile*)fr.At(i); + if (curFile) { + printf("File No. %d: \"%s\"\n", i, curFile->GetName()); + } +// else { +// printf("File No. %d was not loaded. Maximal number of files is %d\n", i, NOCALFILES); +// } + } + + return; +} + +void CsICalib::PrintCuts() { + +// printf("CsICalib::PrintCuts: works wrong\n"); +// 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 + + for (Int_t i = 0; i < cutsCol.GetEntries(); i++) { + curCut = (TCutG*)cutsCol.At(i); + if (curCut) { + printf("Cut No. %d; Name: \"%s\"\n", i, curCut->GetName()); + } /*else { + printf("TOF cut No. %d was not loaded. Maximal number of cuts is %d\n", i, NOCALFILES); + }*/ + } + +/* 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("CsICalib::PrintCuts: End of function.\n"); +*/ + return; +} + +void CsICalib::DrawVariable(const char* variable, Int_t tree, TCanvas *canvas, Int_t lowRange, Int_t upRange) { + +// if (!canvas) TCanvas *c = new TCanvas(); + if (!canvas) return; + + canvas->Clear(); + canvas->Divide(4,4); + + TString canvasTitle; + TString var; + TString con; + + TTree *curTree = 0; + curTree = (TTree*)tr.At(tree); + if (!curTree) { + printf("CsICalib::DrawVariable: Tree No. %d was not found.\n", tree); + return; + } + + canvasTitle.Form("variable: %s; tree: %d", variable, tree); + canvas->SetTitle(canvasTitle.Data()); + + for (Int_t i = 0; i < 16; i++) { + var.Form("%s[%d]", variable, i); + con.Form("%s[%d]>%d && %s[%d]<%d", variable, i, lowRange, variable, i, upRange); + canvas->cd(i+1); + curTree->Draw(var.Data(), con.Data()); + canvas->Update(); + } + +} + +void CsICalib::DrawBeam(TCanvas *canvas, Int_t files, const char* variable) { + + canvas->SetTitle("Beam"); + canvas->Clear(); +// canvas->Divide(files, 2); + canvas->Divide(files, 3); + + TTree *curTree = 0; + TCutG *curCutG = 0; + TString var; + TString con; + + for (Int_t i = 0; i < files; i++) { + canvas->cd(i+1); + curTree = (TTree*)tr.At(i); + if (!curTree) { + printf("CsICalib::DrawBeam: Tree No. %d was not found.\n", i); + continue; + } + curTree->Draw("QDC[0]:TDC[0]", "TDC[0]<1000 && QDC[0]<2000", "cont"); +// 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 < energyPoints; j++) { + if ( cutsCol.At(j) ) { + curCutG = (TCutG*)cutsCol.At(j); + curCutG->Draw("same"); +// printf("CsICalib::DrawBeam: cTOF cut No. %d cannot be drawn, need to repair this function.\n", j); + } + } + + canvas->cd(files+1+i); + curTree->Draw("QDC[0]:QDC[1]", "", "cont"); + for (Int_t j = 0; j < energyPoints; j++) { + if ( cutsCol.At(files + j) ) { + curCutG = (TCutG*)cutsCol.At(files + j); + curCutG->Draw("same"); +// printf("CsICalib::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j); + } + } + + canvas->cd(2*files+1+i); + var.Form("%s[5]:TDC[0]", 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"); +// printf("CsICalib::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j); +// } + } + canvas->Update(); + } +} + +void CsICalib::DrawVariableCut(const char* variable, Int_t tree, TCanvas *canvas, const char* cut1, const char* cut2, Int_t lowRange) { + +// cout << "kjashbdfjka ajsdbf jakhsdb askjdhb" << endl; + + if (!cutsCol.FindObject(cut1)) { + printf("Cut %s was not found.\n", cut1); + return; + } + + canvas->Clear(); + canvas->Divide(4,4); + + TString canvasTitle; + TString var; + TString con; + TString sVariable = variable; + sVariable.ToLower(); + Int_t channel = 0; + +// cout << cutsCol.FindObject(cut1)->GetName() << endl; + + canvasTitle.Form("variable: %s; cut1: %s; cut2: %s; tree: %d", variable, cutsCol.FindObject(cut1)->GetName(), cut2 , tree); + canvas->SetTitle(canvasTitle.Data()); + + TTree *curTree = 0; + curTree = (TTree*)tr.At(tree); + if (!curTree) { + printf("CsICalib::DrawVariableCut: Tree No. %d was not found.\n", tree); + return; + } + + for (Int_t i = 0; i<16; i++) { + if (sVariable.Contains("anc")) { + channel = i+1; + } else { + channel = i; + } +// cout << channel << endl; + var.Form("%s[%d]", variable, channel); + con.Form("%s[%d]>%d", variable, channel, lowRange); + + canvas->cd(i+1); + curTree->SetLineColor(1); + curTree->Draw(var.Data(), con.Data()); + curTree->SetLineColor(3); + con.Form("%s[%d]>0 && %s", variable, channel, cutsCol.FindObject(cut1)->GetName()); + curTree->Draw(var.Data(), con.Data(), "same"); + + TString scut2 = cut2; + if (scut2.Length() != 0) { + if (!cutsCol.FindObject(cut2)) { + printf("Cut %s was not found.\n", cut2); + continue; + } + curTree->SetLineColor(4); + con.Form("%s[%d]>0 && %s", variable, channel, cutsCol.FindObject(cut2)->GetName()); + curTree->Draw(var.Data(), con.Data(), "same"); + } + canvas->Update(); + } + +} + +void CsICalib::GetPeakMean(const char* variable, Int_t tree, Int_t energy, TCanvas *canvas, const char* beamcut, const Int_t nbins, Int_t lowRange) { + + canvas->Clear(); + canvas->Divide(4,4); + + TString var; + TString con; + TString hname; + TString canvasTitle; + + canvasTitle.Form("variable: %s; tree: %d; cut: %s;", variable, tree, cutsCol.FindObject(beamcut)->GetName()); + canvas->SetTitle(canvasTitle.Data()); + + TTree *curTree = 0; + curTree = (TTree*)tr.At(tree); + if (!curTree) { + printf("CsICalib::GetPeakMean: Tree No. %d was not found.\n", tree); + return; + } + + TString sVariable = variable; + sVariable.ToLower(); + Int_t channel = 0; + + for (Int_t i = 0; i<16; i++) { + + if (sVariable.Contains("anc")) { + channel = i+1; + } else { + channel = i; + } + + var.Form("%s[%d]>>hfull[%d][%d]", variable, channel, tree, i); + 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); + 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]->SetLineColor(3); + curTree->Draw(var.Data(), con.Data(), "same"); + + gPad->Update(); + mean[tree][i] = hcut[tree][i]->GetMean(); + meanRMS[tree][i] = hcut[tree][i]->GetRMS(); +// cout << meanRMS[tree][i] << endl; +// cout << hcut[tree][i]->GetMean() << "\t" << hcut[tree][i]->GetRMS() << endl << endl; + + + } + + canvas->Update(); + +} + +void CsICalib::Calibrate(TCanvas *canvas, Bool_t savefile, const char* filename, const char* option) { + + canvas->Clear(); + canvas->Divide(4,4); + +// cout << alphas2.GetSize()+1 << endl; + cout << energyPoints+1 << endl; + + TString gName; + TString gTitle; + +// if (savefile) fGraphs->Open(filename, "RECREATE"); + + TF1 *fnc; + + for (Int_t i = 0; i<16; i++) { + canvas->cd(i+1); + gCal[i] = new TGraphErrors(energyPoints+1); +// FillGraph(gCal[i], energies.GetSize()+1, energies.GetArray(), i); + FillGraph(gCal[i], energyPoints+1, fE.GetArray(), i); +// if (savefile) gCal[i]->Write(); +// gCal[i]->Draw("Al*"); + gCal[i]->Draw("A*"); + gName.Form("g%s%s%d\n", detName.Data(), partName.Data(), i); + gTitle.Form("%s %s\n", detName.Data(), partName.Data()); +// gCal[i]->SetTitle(gTitle.Data()); + gCal[i]->SetName(gName.Data()); + gCal[i]->Fit("pol1"); + fnc = gCal[i]->GetFunction("pol1"); + fnc->SetLineColor(kRed); + fA[i] = fnc->GetParameter(1); + fB[i] = fnc->GetParameter(0); + canvas->Update(); + } + if (savefile) SaveClbGraphs(filename, option); + +} + +void CsICalib::FillGraph(TGraphErrors *g, Int_t npoints, Double_t *energies, Int_t graphNumber, const char* option) { + + TString opt = option; + + //all available energy points and (0,0) + g->SetPoint(0, 0., 0.); + 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]); + g->SetPoint(i+1, mean[i][graphNumber], energies[i]); + g->SetPointError(i+1, meanRMS[i][graphNumber], 0); + } + +// for (Int_t j = 1; j < 4; j++) { +// g->SetPoint(j, 0., 0.); +// } + + +} + +void CsICalib::WriteClbParameters(const char* filename) { + + std::ofstream outfile(filename); + if ( !outfile.is_open() ) { + printf("CsICalib::WriteClbParameters: File %s was not open.\n", filename); + return; + } + + outfile << "#detector:\t" << detName << ",\tparticle:\t" << partName << endl; + outfile << "#channel\tfA\tfB" << endl; + for (Int_t i = 0; i < 16; i++) { + outfile << i << "\t" << fA[i] << "\t" << fB[i] << endl; + } +} + +void CsICalib::SaveClbGraphs(const char* filename, const char* option) { + + cout << "asdasd" << endl; + cout << fGraphs << endl; + if (fGraphs) fGraphs->Close(); + cout << "asdasd" << endl; + fGraphs = new TFile(filename, option); + cout << fGraphs->IsOpen() << endl; + cout << fGraphs->GetName() << endl; + fGraphs->Print(); +// if (!fGraphs->IsOpen()) { +// printf("CsICalib::SaveClbGraphs: file %s was not open.\n", filename); +// return; +// } + + for (Int_t i = 0; i<16; i++) { + fGraphs->cd(); + gCal[i]->Write(); + } + fGraphs->Close(); + return; +} + +void CsICalib::ReadClbParameters(const char* filename) { + + std::ifstream infile(filename); + if ( !infile.is_open() ) { + printf("CsICalib::ReadClbParameters: File %s was not open.\n", filename); + return; + } + + TString line; + Int_t i; + Double_t a, b; + + while (!infile.eof()) { + line.ReadLine(infile); + + if ( line.BeginsWith("#") ) continue; + +// cout << line.Data() << endl; + sscanf(line.Data(), "%d %lf %lf", &i, &a, &b); + fA[i] = a; + fB[i] = b; +// printf("fA[%d]: %f,\tfB[%d]: %f\n", i, fA[i], i, fB[i]); +// printf("fA[%d]: %f,\tfB[%d]: %f\n\n", i, a, i, b); + } + +} + +void CsICalib::DrawClbGraphs(const char* filename, const char* graphname, TCanvas *canvas) { + + printf("CsICalib::DrawClbGraphs: does not work\n"); + return; + + TFile gfile(filename); + + TString gName; +// TGraph *gr; + + for (Int_t i = 0; i < 16; i++) { + gName.Form("%s%d", graphname, i); + gCal[i] = (TGraphErrors*)gfile.Get(gName.Data()); + canvas->cd(i+1); + gCal[i]->Draw("A*"); + } +} + +void CsICalib::DrawEnergyDeposite(const char* variable, TCanvas *canvas, Int_t tree, const char* option) { + if (!canvas) return; + + canvas->Divide(4,4); + + TString opt; + opt = option; + opt.ToLower(); + + TString canvasTitle; + TString var; + TString con; + + TTree *curTree = 0; + curTree = (TTree*)tr.At(tree); + if (!curTree) { + printf("CsICalib::DrawVariable: Tree No. %d was not found.\n", tree); + return; + } + + canvasTitle.Form("variable: %s [MeV]; tree: %d", variable, tree); + canvas->SetTitle(canvasTitle.Data()); + + if (!opt.Contains("same")) canvas->Divide(4,4); + for (Int_t i = 0; i < 16; i++) { + var.Form("%s[%d]*%f + %f", variable, i, fB[i], fA[i]); +// var.Form("%s[%d]*%f + %f", variable, i, fA[i], fB[i]); + con.Form("%s[%d]>0", variable, i); + + if (opt.Contains("same")) { + if (i==0) curTree->Draw(var.Data(), con.Data()); + curTree->Draw(var.Data(), con.Data(), "same"); + } + else { + canvas->cd(i+1); + curTree->Draw(var.Data(), con.Data()); + } + canvas->Update(); + } + + +} diff --git a/Utilities/CsICalib.h b/Utilities/CsICalib.h new file mode 100644 index 0000000000000000000000000000000000000000..42fa9230e0ed61a1a4f36457f04dc56a313d8d8a --- /dev/null +++ b/Utilities/CsICalib.h @@ -0,0 +1,98 @@ +#pragma once + +//#include "TObject.h" +//#include "TROOT.h" +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TCutG.h" +#include "TCanvas.h" +#include "TClonesArray.h" +#include "TH1I.h" +#include "TGraphErrors.h" +#include "TArrayD.h" +#include "TF1.h" + +#define NOCALFILES 5 + +using std::cout; +using std::endl; + +class CsICalib { + private: + + TString detName; + TString partName; + + TClonesArray fr; //TFile +// TFile *fData; + TClonesArray tr; + TClonesArray colFiles; + TClonesArray colTrees; + Int_t nofiles; + TString fileNames[100]; + TString cutNames[100]; + + Int_t energyPoints; + TArrayD fE; + + TFile *fCuts; + TString cutsFileName; + Int_t nCuts; //number of cuts + TClonesArray cutsCol; + + TH1I *hfull[NOCALFILES][16]; //! + TH1I *hcut[NOCALFILES][16]; //! + + Int_t peakMin[NOCALFILES][16]; + Int_t peakMax[NOCALFILES][16]; + Double_t mean[NOCALFILES][16]; + Double_t meanRMS[NOCALFILES][16]; + + TGraphErrors *gCal[16]; + TFile *fGraphs; + TArrayD fA; + TArrayD fB; + + public: +// CsICalib() : a(0), b(0), c(0), p(0){}; + CsICalib(); + CsICalib(const char* parfile); + virtual ~CsICalib(); + + + + void OpenTrees(); + void LoadCuts(); + + void ReadParFile(const char* parfile); + void PrintParameters(const char* option = ""); + void PrintPeakRanges(); + + void DrawVariable(const char* variable, Int_t tree, TCanvas *canvas, Int_t lowRange = 0, Int_t upRange = 4096); + void DrawBeam(TCanvas *canvas, Int_t files, const char* variable); +// void DrawdEE(const char* variable, Int_t tree, TCanvas *canvas); + void DrawVariableCut(const char* variable, Int_t tree, TCanvas *canvas, const char* cut1, const char* cut2 = "", Int_t lowRange = 0); + void GetPeakMean(const char* variable, Int_t tree, 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 FillGraph(TGraphErrors *g, Int_t npoints, Double_t *energies, Int_t graphNumber, const char* option = ""); + void WriteClbParameters(const char* filename); + void SaveClbGraphs(const char* filename, const char* option = "READ"); + void ReadClbParameters(const char* filename); + void DrawClbGraphs(const char* filename, const char* graphname, TCanvas *canvas); + void DrawEnergyDeposite(const char* variable, TCanvas *canvas, Int_t tree, const char* option = ""); + + void PrintTrees(); + void PrintFiles(); + void PrintCuts(); + + Double_t GetA(Int_t i) {return fA[i];}; + Double_t GetB(Int_t i) {return fB[i];}; + + + // Define the class for the cint dictionary + ClassDef (CsICalib,1); +}; diff --git a/Utilities/ReturnCodes.h b/Utilities/ReturnCodes.h new file mode 100644 index 0000000000000000000000000000000000000000..9bf1432ce593c96df4def92a529b9b862612e5f0 --- /dev/null +++ b/Utilities/ReturnCodes.h @@ -0,0 +1,17 @@ +#ifndef RETURN_CODES_H +#define RETURN_CODES_H + +//Some useful return codes + + const static int SUCCESS = 0; + const static int IOEXCEPTION = -1; + const static int NOTFOUND = -2; + const static int NULLPOINTER = -3; + const static int UNKNOWN = -4; + const static int FAILURE = -5; + const static int CDEXCEPTION = -6; + const static int EMPTYCONTAINER = -7; + const static int NOTDEFINED = -8; + const static int EXCEPTION = -11; + +#endif diff --git a/Utilities/Utilities.mk b/Utilities/Utilities.mk new file mode 100644 index 0000000000000000000000000000000000000000..cdbdc4a263f92b4111aa343219766d4f4931847f --- /dev/null +++ b/Utilities/Utilities.mk @@ -0,0 +1,33 @@ +################################################################################ +# AculData input with some variables +################################################################################ + +# UTILITIESLIBS := -lCore -lCint -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf +UTILITIESLIBS := -lCore -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMathCore -lGpad -lGraf +#UTILITIESLIBS := -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic + +# Add inputs and outputs from these tool invocations to the build variables +UTILITIES_HEADERS += \ +$(UTILITIES)/CalPars.h \ +$(UTILITIES)/CsICalib.h \ +$(UTILITIES)/ConfigDictionary.h \ +$(UTILITIES)/linkdef.h + +UTILITIESCPP_SRCS += \ +$(UTILITIES)/CalPars.cpp \ +$(UTILITIES)/CsICalib.cpp \ +$(UTILITIES)/ConfigDictionary.cpp \ +$(UTILITIES)/UtilitiesDict.cpp + +UTILITIESOBJS += \ +$(UTILITIES)/CalPars.o \ +$(UTILITIES)/CsICalib.o \ +$(UTILITIES)/ConfigDictionary.o \ +$(UTILITIES)/UtilitiesDict.o + +UTILITIESCPP_DEPS += \ +$(UTILITIES)/CalPars.d \ +$(UTILITIES)/CsICalib.d \ +$(UTILITIES)/ConfigDictionary.d \ +$(UTILITIES)/UtilitiesDict.d + diff --git a/Utilities/linkdef.h b/Utilities/linkdef.h new file mode 100755 index 0000000000000000000000000000000000000000..c5c0855424b06f3cf8229689e1ba72bb570f1b8b --- /dev/null +++ b/Utilities/linkdef.h @@ -0,0 +1,11 @@ +#ifdef __CINT__ +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class CalPars; +#pragma link C++ class CsICalib; +#pragma link C++ class ConfigDictionary; + +#endif + diff --git a/html.cxx b/html.cxx deleted file mode 100755 index acc65a4861d28f0afd8fea01ab663ee1e319d421..0000000000000000000000000000000000000000 --- a/html.cxx +++ /dev/null @@ -1,11 +0,0 @@ -{ - gSystem->Load("./libAculData.so"); - gSystem->Load("./libTELoss.so"); - - THtml h; - - h.SetInputDir("."); - - h.MakeAll(); - -} diff --git a/htmldoc.cxx b/htmldoc.cxx new file mode 100755 index 0000000000000000000000000000000000000000..354a52332b42b2b60df38ea487897af7d844b443 --- /dev/null +++ b/htmldoc.cxx @@ -0,0 +1,18 @@ +#include "THtml.h" + +void htmldoc() +{ + gSystem->Load("./libTELoss.so"); + gSystem->Load("./libAculCalib.so"); + gSystem->Load("./libAculData.so"); + + THtml h; + + h.SetInputDir("."); + h.SetOutputDir("htmldoc"); + h.SetProductName("AculUtils"); + + + h.MakeAll(); + +} diff --git a/macros/calibration_CsI/calibrationSQ13Alpha.cxx b/macros/calibration_CsI/calibrationSQ13Alpha.cxx index d76898a24b3595e9055b0d56dca9305fba277a56..972b97f5e3bbef5bbe2db573a4d867a0839db607 100644 --- a/macros/calibration_CsI/calibrationSQ13Alpha.cxx +++ b/macros/calibration_CsI/calibrationSQ13Alpha.cxx @@ -33,8 +33,14 @@ // cal.DrawBeam(myCanv, 4, "SQ13"); //return; +// <<<<<<< HEAD // cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16"); // cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp"); +// ======= + cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16"); + return; + cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp"); +// >>>>>>> newArch //return; // cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp"); // cal.DrawVariableCut("SQ13", 2, c3, "cutSQ13Alpha26", "cutSQ13Alpha26Amp"); diff --git a/macros/myMacros/SQ13Alpha.par b/macros/myMacros/SQ13Alpha.par new file mode 100644 index 0000000000000000000000000000000000000000..c11a380d0e2ff7fcef3c24e783868d8eebc44b62 --- /dev/null +++ b/macros/myMacros/SQ13Alpha.par @@ -0,0 +1,103 @@ +energies 4 #number of energy points +crystals 16 #number of crystals in detector +files 4 +#fileID fileName +#fileID not implemented yet +data/csi_13_Ealpha16.root +data/csi_13_Ealpha21.root +data/csi_13_Ealpha26.root +data/csi_13_Ealpha30.root +#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 +3 1200 1400 +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 +13 1140 1380 +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 +7 1450 1700 +8 1150 1400 +9 900 1100 +10 1350 1600 +11 1250 1500 +12 950 1200 +13 1600 1860 +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 +3 2250 2500 +4 1720 2000 +5 1200 1450 +6 1160 1400 +7 1800 2200 +8 1600 1850 +9 1200 1400 +10 1750 2100 +11 1650 1900 +12 1250 1600 +13 2150 2450 +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 +3 2590 3000 +4 2100 2350 +5 1500 1650 +6 1400 1600 +7 2250 2550 +8 1900 2150 +9 1400 1600 +10 2050 2600 +11 1950 2200 +12 1580 1710 +13 2550 2900 +14 1100 1250 +15 1800 2000 \ No newline at end of file diff --git a/macros/myMacros/cuts b/macros/myMacros/cuts new file mode 120000 index 0000000000000000000000000000000000000000..76679909cc6ebded0c3b6c69222f1b126f59fe34 --- /dev/null +++ b/macros/myMacros/cuts @@ -0,0 +1 @@ +../calibration_CsI/cuts/ \ No newline at end of file diff --git a/macros/myMacros/data b/macros/myMacros/data new file mode 120000 index 0000000000000000000000000000000000000000..3080039a9937febeeda9884323fa41947b295ef1 --- /dev/null +++ b/macros/myMacros/data @@ -0,0 +1 @@ +../calibration_CsI/data/ \ No newline at end of file diff --git a/macros/myMacros/parTest.cxx b/macros/myMacros/parTest.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f82aa3a91cca32300720327fa846c85810418fef --- /dev/null +++ b/macros/myMacros/parTest.cxx @@ -0,0 +1,83 @@ +//#include "../../AculCalib/AculCalPars.h" +//#include "../../AculCalib/AculCalibScint.h" +//#include "../../AculCalib/AculCalibSi.h" +//#include "TSystem.h" +//#include +// +using std::cout; +using std::endl; + +void parTest() +{ + gSystem->Load("../../libTELoss.so"); + gSystem->Load("../../libAculCalib.so"); + +/* AculCalibSi ps; + ps.SetParFile("parforcal.par"); + ps.Init(); + ps.PrintParameters(); + return; +*/ + +// AculCalParsScintFile pp; +// pp.SetNoCrystals(16); +// pp.Init(); +// return; + + cout << "-----------------------" << endl; + + AculCalParsScint p; + p.SetParFile("SQ13Alpha.par"); + p.Init(); +// cout << p.GetFileName(0) << endl; +// cout << p.GetFileName(1) << endl; +// 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(""); +// return; + + cout << "-----------------------" << endl; + + AculCalibScint c; + c.SetParFile("SQ13Alpha.par"); + c.Init(); + c.PrintParameters(); +// return; + +// c.PrintParameters("all"); + +// c.OpenFiles(); + +// c.LoadTrees(); + c.PrintTrees(); + + TCanvas *c1 = new TCanvas("c1", "Plain"); + TCanvas *c2 = new TCanvas("c2", "Plain"); + TCanvas *c3 = new TCanvas("c3", "Plain"); + TCanvas *c4 = new TCanvas("c4", "Plain"); +// c.DrawVariable("SQ13", 2, c1); +// c.DrawBeam(c1, 4, "SQ13"); + +// 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(); +// c.PrintCuts(); + +// TCanvas *c1 = new TCanvas("c1", "Plain"); +// c.GetPeakMean("SQ13", 0, 0, c1, "cutSQ13Alpha16Amp", 256); + + TCanvas *cCal = new TCanvas("cCal", "calibration Alpha"); + c.Calibrate(cCal); + + return; +} diff --git a/makefile b/makefile index 92873f21f70ff41970df852ba21f8116cc118009..5d254efd0f5855995fd1ca83061be217a6346768 100755 --- a/makefile +++ b/makefile @@ -25,25 +25,41 @@ ROOTCFLAGS = $(shell root-config --cflags) PWD = $(shell pwd) #INSTALLFOLDER = $(HOME)/AculLib +UTILITIES = $(PWD)/Utilities ACULDATA = $(PWD)/AculData +ACULCALIB = $(PWD)/AculCalib TELOSS = $(PWD)/TELoss +-include $(UTILITIES)/Utilities.mk -include $(ACULDATA)/AculData.mk +-include $(ACULCALIB)/AculCalib.mk -include $(TELOSS)/TELoss.mk -all: libAculData.so \ +all: libUtilities.so \ + libAculData.so \ + libAculCalib.so \ libTELoss.so #ROOT html documentation, it will be done as a program which will be alsa compiled by this makefile, program will be as a last condition after all of the libraries -htmldoc: libAculData.so +htmldoc: libAculData.so libAculCalib.so libTELoss.so -$(RM) htmldoc - root -l -q html.cxx + root -l -q htmldoc.cxx + +clean_htmldoc: + -$(RM) htmldoc + -@echo ' ' clean: + -$(RM) $(UTILITIESOBJS) $(UTILITIESCPP_DEPS) + -$(RM) $(UTILITIES)/UtilitiesDict.* libUtilities.so + -@echo ' ' -$(RM) $(ACULDATAOBJS) $(ACULDATACPP_DEPS) -$(RM) $(ACULDATA)/AculDataDict.* *.pcm $(ACULDATA)/*.pcm -$(RM) libAculData.so -@echo ' ' + -$(RM) $(ACULCALIBOBJS) $(ACULCALIBCPP_DEPS) + -$(RM) $(ACULCALIB)/AculCalibDict.* libAculCalib.so + -@echo ' ' -$(RM) $(TELOSSOBJS) $(TELOSSCPP_DEPS) -$(RM) $(TELOSS)/TELossDict.* *.pcm $(TELOSS)/*.pcm -$(RM) libTELoss.so @@ -51,13 +67,28 @@ clean: -$(RM) htmldoc -@echo ' ' +# <<<<<<< HEAD # Those *Dictionary* files below need special treating: +$(UTILITIES)/UtilitiesDict.cpp: + -@echo 'Pre-building UtilitiesDict.cpp and UtilitiesDict.h files' + -rootcling -f $(UTILITIES)/UtilitiesDict.cpp -p $(UTILITIES_HEADERS) + -@echo 'Creating: link to UtilitiesDict_rdict.pcm' + -ln -s $(UTILITIES)/UtilitiesDict_rdict.pcm . + -@echo ' ' + $(ACULDATA)/AculDataDict.cpp: -@echo 'Pre-building AculDataDict.cpp and AculDataCint.h files' -rootcling -f $(ACULDATA)/AculDataDict.cpp -p $(ACULDATA_HEADERS) -@echo 'Creating: link to AculDataDict_rdict.pcm' -ln -s $(ACULDATA)/AculDataDict_rdict.pcm . -@echo ' ' + +$(ACULCALIB)/AculCalibDict.cpp: + -@echo 'Pre-building AculCalibDict.cpp and AculCalibDict.h files' + -rootcling -f $(ACULCALIB)/AculCalibDict.cpp -p $(ACULCALIB_HEADERS) + -@echo 'Creating: link to AculCalibDict_rdict.pcm' + -ln -s $(ACULDATA)/AculCalibDict_rdict.pcm . + -@echo ' ' $(TELOSS)/TELossDict.cpp: -@echo 'Pre-building TELossDict.cpp and TELossDict_rdict.pcm files' @@ -67,13 +98,49 @@ $(TELOSS)/TELossDict.cpp: -ln -s $(TELOSS)/TELossDict_rdict.pcm . -@echo ' ' +# ======= +# Those *Cint* files below need special treating: +# $(UTILITIES)/UtilitiesDict.cpp: +# -@echo 'Pre-building UtilitiesDict.cpp and UtilitiesDict.h files' +# -rootcint -f $(UTILITIES)/UtilitiesDict.cpp -c -p $(UTILITIES_HEADERS) +# -@echo ' ' + +# $(ACULDATA)/AculDataCint.cpp: +# -@echo 'Pre-building AculDataCint.cpp and AculDataCint.h files' +# -rootcint -f $(ACULDATA)/AculDataCint.cpp -c -p $(ACULDATA_HEADERS) +# -@echo ' ' + +# $(ACULCALIB)/AculCalibDict.cpp: +# -@echo 'Pre-building AculCalibDict.cpp and AculCalibDict.h files' +# -rootcint -f $(ACULCALIB)/AculCalibDict.cpp -c -p $(ACULCALIB_HEADERS) +# -@echo ' ' + +# $(TELOSS)/TELossCint.cpp: +# -@echo 'Pre-building TELossCint.cpp and TELossCint.h files' +# -rootcint -f $(TELOSS)/TELossCint.cpp -c -p $(TELOSS)/TELoss.h $(TELOSS)/linkdef.h +# >>>>>>> newArch + #*.so files +libUtilities.so: $(UTILITIESOBJS) + @echo 'Building target: $@' + @echo 'Invoking: $(CC) Linker' + $(CC) -L $(ROOTLIBS) -shared -o"libUtilities.so" $(UTILITIESOBJS) $(UTILITIESLIBS) + @echo 'Finished building target: $@' + @echo ' ' + libAculData.so: libTELoss.so $(ACULDATAOBJS) @echo 'Building target: $@' @echo 'Invoking: GCC C++ Linker' $(CC) -L . -L $(ROOTLIBS) -shared -o"libAculData.so" $(ACULDATAOBJS) $(ACULDATALIBS) @echo 'Finished building target: $@' @echo ' ' + +libAculCalib.so: libTELoss.so $(ACULCALIBOBJS) + @echo 'Building target: $@' + @echo 'Invoking: GCC C++ Linker' + $(CC) -L . -L $(ROOTLIBS) -shared -o"libAculCalib.so" $(ACULCALIBOBJS) $(ACULCALIBLIBS) + @echo 'Finished building target: $@' + @echo ' ' libTELoss.so: $(TELOSSOBJS) @echo 'Building target: $@'