...
 
Commits (9)
/*
* 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();
}
/*
* AculCalParameters.h
*
* Created on: Oct 20, 2016
* Author: vratik
*/
#ifndef ACULCALIB_ACULCALPARS_H_
#define ACULCALIB_ACULCALPARS_H_
#include <iostream>
#include <fstream>
#include <vector>
#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<Double_t> 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_ */
/*
* 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);
}
}*/
/*
* 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<AculCalParsScintFile> fFilePars;
//todo following parameters should be probably moved to AculCalParsScintFile
//as far as they are related to each file separately
// vector<TString> fCutName;
// vector<TCutG> 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_ */
/*
* AculCalParsScintFile.cpp
*
* Created on: Oct 25, 2016
* Author: vratik
*/
#include "AculCalParsScintFile.h"
AculCalParsScintFile::AculCalParsScintFile() {
// TODO Auto-generated constructor stub
// fFileID = 0;
fNoCrystals = 0;
fE = 0;
Init();
}
AculCalParsScintFile::~AculCalParsScintFile() {
// TODO Auto-generated destructor stub
}
void AculCalParsScintFile::Init() {
// fNoCrystals = 0;
// fE = 0;
fPeakMin.resize(fNoCrystals);
fPeakMax.resize(fNoCrystals);
fCutName.resize(0);
}
const char* AculCalParsScintFile::GetCutName(Int_t i) {
if ( i >= (Int_t)fCutName.size() ) {
cout << "\"AculCalParsScintFile::GetCutName\" Cut number " << i
<< " is not available. Maximal i can be set as "
<< fCutName.size() - 1 << "." << endl;
}
return fCutName[i].Data();
}
TCutG* AculCalParsScintFile::GetCut(Int_t i) {
if ( i >= (Int_t)fCut.size() ) {
cerr << "\"AculCalParsScintFile::GetCut\" index i cannot be higher than " << fCut.size() - 1 << endl;
return 0;
}
return &fCut[i];
}
TCutG* AculCalParsScintFile::GetCut(const char* cutName) {
const TString cName = cutName;
for (Int_t i = 0; i <= (Int_t)fCut.size(); i++) {
if (cName.EqualTo(fCutName[i])) { return &fCut[i]; }
}
cerr << "\"AculCalParsScint::GetCut\" cut \"" << cutName << "\" was not found." << endl;
return 0;
}
Int_t AculCalParsScintFile::GetPeakMin(Int_t crystalID) {
if ( crystalID >= (Int_t)fPeakMin.size() ) {
cerr << "\"AculCalParsScintFile::GetPeakMin\" index \"crystal\" cannot be higher than "
<< fPeakMin.size() - 1 << endl;
return 0;
}
return fPeakMin[crystalID];
}
Int_t AculCalParsScintFile::GetPeakMax(Int_t crystalID) {
if ( crystalID >= (Int_t)fPeakMax.size() ) {
cerr << "\"AculCalParsScintFile::GetPeakMax\" index \"crystal\" cannot be higher than "
<< fPeakMax.size() - 1 << endl;
return 0;
}
return fPeakMax[crystalID];
}
void AculCalParsScintFile::SetRawFileName(const char* filename) {
fFileName = filename;
// fFileID++;
}
void AculCalParsScintFile::SetCutFileName(const char* filename) {
fCutsFileName = filename;
}
void AculCalParsScintFile::SetNoCuts(Int_t noCuts) {
fCutName.resize(noCuts);
}
void AculCalParsScintFile::SetCutName(Int_t i, const char* cutname) {
if ( i >= (Int_t)fCutName.size() ) {
cout << "\"AculCalParsScintFile::SetCutName\" vector for cut names is not enough large. Maximal i can be set as "
<< fCutName.size() - 1 << "." << endl;
}
fCutName[i] = cutname;
}
void AculCalParsScintFile::SetPeakMin(Int_t crystalID, Double_t limit) {
if ( crystalID>=(Int_t)fPeakMin.size() ) {
cerr << "\"AculCalParsScintFile::SetPeakMin\" maximal allowed crystal ID is " << fPeakMin.size() << "."<< endl;
return;
}
fPeakMin[crystalID] = limit;
}
void AculCalParsScintFile::SetPeakMax(Int_t crystalID, Double_t limit) {
if ( crystalID>=(Int_t)fPeakMax.size() ) {
cerr << "\"AculCalParsScintFile::SetPeakMax\" maximal allowed crystal ID is " << fPeakMax.size() << "."<< endl;
return;
}
fPeakMax[crystalID] = limit;
}
void AculCalParsScintFile::LoadCuts() {
if (fCutsFileName.Length() == 0) {
printf("AculCalParsScintFile::LoadCuts: Name of file (*.root) with cuts was not provided.\n");
printf("AculCalParsScintFile::LoadCuts: No cuts has been loaded.\n");
return;
}
TFile cutFile(fCutsFileName.Data(), "READ");
if(!cutFile.IsOpen()) {
cerr << "\"AculCalParsScintFile::LoadCuts\" File " << fCutsFileName.Data()
<< " was not open and no cuts were loaded." << endl;
return;
}
for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) {
TCutG *currentCut = (TCutG*)cutFile.Get(fCutName[i]);
if (!currentCut) {
cout << "\"AculCalParsScintFile::LoadCuts\" Cut \"" << fCutName[i]
<< "\" was not found in file " << fCutsFileName << "." << endl;
continue;
}
fCut.push_back(*currentCut);
}
}
/*void AculCalParsScintFile::Reset() {
}*/
/*
* AculCalParsScintFile.h
*
* Created on: Oct 25, 2016
* Author: vratik
*/
#ifndef ACULCALIB_ACULCALPARSSCINTFILE_H_
#define ACULCALIB_ACULCALPARSSCINTFILE_H_
#include <vector>
#include <iostream>
#include "TFile.h"
#include "TCutG.h"
#include "TString.h"
using std::vector;
using std::cout;
using std::endl;
using std::cerr;
class AculCalParsScintFile {
private:
// UInt_t fFileID;
Int_t fNoCrystals;
Double_t fE;
// TCutG fCuts;
TString fFileName;
TString fCutsFileName;
vector<TString> fCutName;
vector<TCutG> fCut;
vector<Int_t> fPeakMin;
vector<Int_t> fPeakMax;
public:
AculCalParsScintFile();
virtual ~AculCalParsScintFile();
ClassDef(AculCalParsScintFile, 1);
void Init();
Int_t GetNoCrystals() {return fNoCrystals;}
const char* GetRawFileName() {return fFileName.Data();}
const char* GetCutFileName() {return fCutsFileName.Data();}
Double_t GetEnergy() {return fE;}
Int_t GetNoCuts() {return (Int_t)fCut.size();}
const char* GetCutName(Int_t i);
TCutG* GetCut(Int_t i);
TCutG* GetCut(const char* cutName);
Int_t GetPeakMin(Int_t crystalID);
Int_t GetPeakMax(Int_t crystalID);
void SetNoCrystals(Int_t nocrystals) {fNoCrystals = nocrystals;};
void SetEnergy(Double_t energy) {fE = energy;};
void SetRawFileName(const char* filename);
void SetCutFileName(const char* filename);
void SetNoCuts(Int_t noCuts);
void SetCutName(Int_t i, const char* cutname);
void SetPeakMin(Int_t crystalID, Double_t limit);
void SetPeakMax(Int_t crystalID, Double_t limit);
void LoadCuts();
// void Reset();
private:
// void LoadCuts();
};
#endif /* ACULCALIB_ACULCALPARSSCINTFILE_H_ */
/*
* AculCalParsSi.cpp
*
* Created on: Oct 26, 2016
* Author: vratik
*/
#include "AculCalParsSi.h"
AculCalParsSi::AculCalParsSi() {
// TODO Auto-generated constructor stub
Reset();
}
AculCalParsSi::AculCalParsSi(const char* parFile) {
SetParFile(parFile);
Init();
}
AculCalParsSi::~AculCalParsSi() {
// TODO Auto-generated destructor stub
}
void AculCalParsSi::Init() {
SetELosses();
SetPars();
SetCalEnergies();
}
void AculCalParsSi::Reset() {
fParFileName = "";
kRaNOPEAKS = 0;
fEnergyInput.clear();
fLowerChannel = 0.;
fUpperChannel = 0.;
fLowerPeakRelativeHight = 0.;
fUpperPeakRelativeHight = 0.;
fPeakPositionTolerance = 0.;
fFitFuncLineWidth = 1;
fFitMinSigma = 0.;
fFitPeakThreshold = 0.;
fDeadLayer = 0.;
return;
}
void AculCalParsSi::SetPars() {
if (fParFileName.Length()==0) {
cerr << "\'AculCalibration::SetInputsParameters\" File with input parameters was not set." << endl;
return;
}
const Int_t lineLength = 400;
Char_t line[lineLength];
Char_t parameter[100];
Char_t identificator[100];
ifstream fipr;
fipr.open(fParFileName.Data());
if (!fipr.is_open()) {
cerr << "\"AculCalibration::SetInputsParameters\" File with input parameters \""
<< fParFileName << "\" was not opened." << endl;
return;
}
cout << "\"AculCalibration::SetInputsParameters\" File with input parameters \""
<< fParFileName << "\" will be processed." << endl;
while (!fipr.eof()) {
fipr.getline(line, lineLength);
if (strlen(line) < 2) {
continue;
}
sscanf(line, "%s %s", parameter, identificator);
if ( strcmp(identificator, "noPeaks") == 0 ) {
kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
fEnergyInput.resize(kRaNOPEAKS);
fEtab.resize(kRaNOPEAKS);
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
fipr.getline(line, lineLength);
sscanf(line, "%s", parameter);
fEnergyInput[i] = static_cast<Double_t>(atof(parameter));
}
continue;
}//if
if ( strcmp(identificator, "lowerChannel") == 0 ) {
sscanf(line, "%s", parameter);
fLowerChannel = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "upperChannel") == 0 ) {
sscanf(line, "%s", parameter);
fUpperChannel = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "lowerPeakHight") == 0 ) {
sscanf(line, "%s", parameter);
fLowerPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "upperPeakHight") == 0 ) {
sscanf(line, "%s", parameter);
fUpperPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "peakPositionTolerance") == 0 ) {
sscanf(line, "%s", parameter);
fPeakPositionTolerance = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fitFunctionLineWidth") == 0 ) {
sscanf(line, "%s", parameter);
fFitFuncLineWidth = static_cast<Width_t>(atoi(parameter));
}
if ( strcmp(identificator, "minFitSigma") == 0 ) {
sscanf(line, "%s", parameter);
fFitMinSigma = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fitHightThreshold") == 0 ) {
sscanf(line, "%s", parameter);
fFitPeakThreshold = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "deadLayer") == 0 ) {
sscanf(line, "%s", parameter);
fDeadLayer = static_cast<Double_t>(atof(parameter));
}
}
fipr.close();
return;
}
void AculCalParsSi::PrintParameters(const char* option)
{
//print alpha source parameters
cout << "AculCalibration::PrintInputParameters:" << endl;
cout << "\tNumber of peaks: " << kRaNOPEAKS << endl;
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
cout << "\t\tfEnergyInput[" << i << "] = " << fEnergyInput[i] << " MeV" << endl;
}
cout << "\tEnergies used for calibration:" << endl;
cout << "\t(deadLayer: " << fDeadLayer << " mcm)" << endl;
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
cout << "\t\tfE[" << i << "] = " << fEtab[i] << " MeV" << endl;
}
cout << "\tlowerChannel: " << fLowerChannel << "; upperChannel: " << fUpperChannel << ";" << endl;
cout << "\tlowerPeakHight: " << fLowerPeakRelativeHight << "; upperPeakHight: " << fUpperPeakRelativeHight << ";" << endl;
cout << "\tfitHightThreshold: " << fFitPeakThreshold << "; minFitSigma: " << fFitMinSigma << ";" << endl;
cout << "\tpeakPositionTolerance: " << fPeakPositionTolerance << ";" << endl;
cout << "\tfitFunctionLineWidth: " << fFitFuncLineWidth << ";" << endl;
return;
}
void AculCalParsSi::SetELosses() {
Info("AculCalibration::SetELosses", "Combination of aplha particle with silicon material only.");
fAlphaSi.SetEL(1, 2.321); // density in g/cm3
fAlphaSi.AddEL(14., 28.086, 1); //Z, mass
// mSi.SetZP(1., 1.); //protons
fAlphaSi.SetZP(2., 4.); //alphas, Z, A
fAlphaSi.SetEtab(100000, 200.); // ?, MeV calculate ranges
fAlphaSi.SetDeltaEtab(300);
}
void AculCalParsSi::SetCalEnergies() {
if (fDeadLayer<=0.) {
Warning("AculCalibration::SetCalEnergies", "Dead layer was set equal or less than 0.");
for(Int_t i = 0; i < kRaNOPEAKS; i++) {
fEtab[i] = fEnergyInput[i];
}
Info("AculCalibration::SetCalEnergies", "Energies used for calibration are the same as input file.");
return;
}
for(Int_t i = 0; i < kRaNOPEAKS; i++) {
fEtab[i] = fAlphaSi.GetE(fEnergyInput[i], fDeadLayer);
}
Info("AculCalibration::SetCalEnergies", "Energies used for calibration considering %f mcm dead layer were set.", fDeadLayer);
return;
}
/*
* AculCalParsSi.h
*
* Created on: Oct 26, 2016
* Author: vratik
*/
#ifndef ACULCALIB_ACULCALPARSSI_H_
#define ACULCALIB_ACULCALPARSSI_H_
#include "AculCalPars.h"
#include "../TELoss/TELoss.h"
class AculCalParsSi: public AculCalPars {
private:
Int_t kRaNOPEAKS;
// TArrayD fEnergy; //energy after passing through deadlayer
vector<Double_t> fEnergyInput; //incidental energy, set from .par
Double_t fLowerChannel;
Double_t fUpperChannel;
Double_t fLowerPeakRelativeHight; //pouziva se, private
Double_t fUpperPeakRelativeHight; //pouziva se, private, nastavit nenulovou prednastavenou hodnotu
Double_t fPeakPositionTolerance; //pouziva se, private
Width_t fFitFuncLineWidth; //private
Double_t fFitMinSigma; //pouziva se, private
Double_t fFitPeakThreshold; //pouziva se, private, prozkoumat, k cemu vlastne slouzi ve fci ShowPeaks, popremyslet o vhodnem prednastaveni v konstruktoru
Double_t fDeadLayer; //in mcm; it will be used for recalculation of energies outcomming from the alpha source
TELoss fAlphaSi;
public:
AculCalParsSi();
AculCalParsSi(const char* parFile);
virtual ~AculCalParsSi();
ClassDef(AculCalParsSi, 1);
void Init();
// Function which loads text file containing parameters for calibration
//
// -inputparfile: file containing information on calibration source
// -In file with the data must be preserved systematic
// -There cannot be whitelines (probably)
//
//Example of "parforcal.par" file
//................................................................
//223Ra
//
//4 nopeaks //number of peaks
//4.415 E1 //in MeV
//5.153 E2 //in MeV
//5.683 E3 //in MeV
//7.419 E4 //in MeV
//100 lowerchannel //in channels
//4096 upperchannel //in channels
//0.1 lowerpeakhight //in relative units
//0.1 upperpeakhight //in relative units
//0.1 peakpositiontolerance //in relative units
//2 fitfunctionlinewidth //integer 1 - 10
//5 minfitsigma //minimal sigma of the peaks to be fitted
//0.4 fithightthreshold //
//................................................................
void Reset();
void PrintParameters(const char* option = "");
//I hope this is selfunderstanding function
private:
void SetPars();
void SetELosses();
void SetCalEnergies();
};
#endif /* ACULCALIB_ACULCALPARSSI_H_ */
#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);
////////////////////////////////////////////////////////////////
// //
// TMyClass //
// //
// This is the description block. //
// //
////////////////////////////////////////////////////////////////
#ifndef ACULCALIB_ACULCALIB_H_
#define ACULCALIB_ACULCALIB_H_
//#include "TObject.h"
//#include "TROOT.h"
#include <iostream>
#include <fstream>
#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<Double_t> fA; //calibration parameter A
vector<Double_t> 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
################################################################################
# AculData input with some variables
################################################################################
ACULCALIBLIBS := -lCore -lCint -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)/AculCalibCint.cpp
ACULCALIBOBJS += \
$(ACULCALIB)/AculCalib.o \
$(ACULCALIB)/AculCalibScint.o \
$(ACULCALIB)/AculCalibSi.o \
$(ACULCALIB)/AculCalPars.o \
$(ACULCALIB)/AculCalParsScint.o \
$(ACULCALIB)/AculCalParsScintFile.o \
$(ACULCALIB)/AculCalParsSi.o \
$(ACULCALIB)/AculCalibCint.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)/AculCalibCint.d
\ No newline at end of file
#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();
}
//////////////////////////////////////////////////////////////////////////////
// //
// 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 <iostream>
#include <fstream>
//#include <vector>
#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<TH1I*> > fHistFull;
vector< vector<TH1I*> > fHistCut;
//
vector< vector<Double_t> > fMeanPeakPos; //doc of class member
vector< vector<Double_t> > 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();
};
/*
* 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());
}
/*
* 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_ */
#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
#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 <nopars; i++) {
line.ReadLine(infile);
// if ( line.BeginsWith("#") ) continue;
// cout << line.Data() << endl;
sscanf(line.Data(), "%lf %lf %d", &b, &a, &c);
fA[i] = a;
fB[i] = b;
fC[i] = c;
// 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 CalPars::PrintClbPars() {
printf("CalPars::PrintClbPars: %d calibration parameters from %s file.\n", fA.GetSize(), clbFile.Data());
for (Int_t i = 0; i<fA.GetSize(); i++) {
printf("%d\t%f\t%f\t%f\n", i, GetA(i), GetB(i), GetC(i));
}
}
#pragma once
#include <iostream>
#include <fstream>
#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);
};
#include "ConfigDictionary.h"
ClassImp(ConfigDictionary);
//////////////////////////////////////////////////////////////////////////////
// BEGIN_HTML
// <p><font size="4"><b>Config(uration)Dictionary class</b></font></p>
// <br>
// <i>Author: Bartlomiej Hnatio 2012-08-06</i>
// <br><br>
// 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.
// <br><br><br>
// <u>Most simple example of usage:</u>
// <br><br>Suppose you have two variables fD and fI:
// <pre>-------------------
// Double_t fD;
// Int_t fI;
//
//fD = 3.14;
//fI = 2012;
//-----------------------</pre>
// To save its parameters into file you can create ConfigDictionary class
// instance and use <b>SetDouble</b> and <b>SetInt</b> functions to
// insert parameters values with arbitrarly defined keys (let them be
// really long and descriptive in this example):
// <pre>---------------------
//ConfigDictionary CD;
//CD.SetDouble("Most important variable D",fD);
//CD.SetInt("Equally important integer I",fI);
//---------------------</pre>
// Now configuration string saved in <b>ConfigDictionary</b> CD
// can be obtained with <b>ToString</b> method and should look like:
//<pre>
//"Most important variable D"="3.14" "Equally important integer I"="2012"
//</pre>
// <br> It can be easily saved to a file using simplest <i>fstream</i>
// 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.
// <br>
// <br>
// 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:
// <pre>-----------------
//"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"
//(...)
//----------------</pre>
// 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 <b>FromString</b> method. After this, you can
// extract parameters (using <b>GetDouble</b> and <b>GetInt</b> methods)
// with using same keys as were used for saving, and ConfigDictionary
// will return their values:
// <pre>----------------------
//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");
//--------------------</pre>
//<br>
// 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 <b>Get</b> 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 <b>Get</b> 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:
//<pre>-------------------
//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!");
//}
//--------------------</pre>
//
// 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<std::string,std::string>::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<<"\""<<it->first<<"\""<<"="<<"\""<<it->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<int>(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<double>(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<double>(key,value);
}
//_____________________________________________________________________________
void ConfigDictionary::SetInt(std::string key,int value){
//Sets value to key, converts int to string first
Set<int>(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";
}
#ifndef CONFIGDICTIONARY_H
#define CONFIGDICTIONARY_H
#include "TObject.h"
#include "ReturnCodes.h"
#include "TStopwatch.h"
#include "TROOT.h"
#include <string>
#include <map>
#include <sstream>
#include <vector>
#include "TLorentzVector.h"
#include <iostream>
#include "TMath.h"
#include "TString.h"
#include "TTree.h"
#define ARRAY_DELIM '$'
class ConfigDictionary{
public:
typedef std::map<std::string,std::string>::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 <class T>
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 <class T>
void Set(std::string key,T value){
//Sets value to key, converts T to string first with << operator.
std::stringstream ss;
ss<<value;
configMap[key] = ss.str();
}
//_________________________________________________________________________
template <class T>
int GetArray(std::string key,std::vector<T> & 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 <class T>
void SetArray(std::string key,std::vector<T> & input){
std::ostringstream arrayString;
for (unsigned int i=0;i<input.size();i++){
arrayString<<input[i]<<ARRAY_DELIM;
}
configMap[key] = arrayString.str();
}
private:
std::map<std::string,std::string> configMap;
};
#endif
#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();
}
}
#pragma once
//#include "TObject.h"
//#include "TROOT.h"
#include <iostream>
#include <fstream>
#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);
};
#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
################################################################################
# AculData input with some variables
################################################################################
UTILITIESLIBS := -lCore -lCint -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)/UtilitiesCint.cpp
UTILITIESOBJS += \
$(UTILITIES)/CalPars.o \
$(UTILITIES)/CsICalib.o \
$(UTILITIES)/ConfigDictionary.o \
$(UTILITIES)/UtilitiesCint.o
UTILITIESCPP_DEPS += \
$(UTILITIES)/CalPars.d \
$(UTILITIES)/CsICalib.d \
$(UTILITIES)/ConfigDictionary.d \
$(UTILITIES)/UtilitiesCint.d
#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
#include "THtml.h"
void htmldoc()
{ {
gSystem->Load("./libAculData.so");
gSystem->Load("./libTELoss.so"); gSystem->Load("./libTELoss.so");
gSystem->Load("./libAculCalib.so");
gSystem->Load("./libAculData.so");
THtml h; THtml h;
h.SetInputDir("."); h.SetInputDir(".");
h.SetOutputDir("htmldoc");
h.SetProductName("AculUtils");
h.MakeAll(); h.MakeAll();
......
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
// cal.DrawBeam(c4, 4, "SQ13"); // cal.DrawBeam(c4, 4, "SQ13");
//return; //return;
// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16"); // cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16");
// return;
// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp"); // cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp");
//return; //return;
// cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp"); // cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp");
......
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
../calibration_CsI/cuts/
\ No newline at end of file
../calibration_CsI/data/
\ No newline at end of file
//#include "../../AculCalib/AculCalPars.h"
//#include "../../AculCalib/AculCalibScint.h"
//#include "../../AculCalib/AculCalibSi.h"
//#include "TSystem.h"
//#include <iostream>
//
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;
}
...@@ -25,24 +25,40 @@ ROOTCFLAGS = $(shell root-config --cflags) ...@@ -25,24 +25,40 @@ ROOTCFLAGS = $(shell root-config --cflags)
PWD = $(shell pwd) PWD = $(shell pwd)
#INSTALLFOLDER = $(HOME)/AculLib #INSTALLFOLDER = $(HOME)/AculLib
UTILITIES = $(PWD)/Utilities
ACULDATA = $(PWD)/AculData ACULDATA = $(PWD)/AculData
ACULCALIB = $(PWD)/AculCalib
TELOSS = $(PWD)/TELoss TELOSS = $(PWD)/TELoss
-include $(UTILITIES)/Utilities.mk
-include $(ACULDATA)/AculData.mk -include $(ACULDATA)/AculData.mk
-include $(ACULCALIB)/AculCalib.mk
-include $(TELOSS)/TELoss.mk -include $(TELOSS)/TELoss.mk
all: libAculData.so \ all: libUtilities.so \
libAculData.so \
libAculCalib.so \
libTELoss.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 #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 -$(RM) htmldoc
root -l -q html.cxx root -l -q htmldoc.cxx
clean_htmldoc:
-$(RM) htmldoc
-@echo ' '
clean: clean:
-$(RM) $(UTILITIESOBJS) $(UTILITIESCPP_DEPS)
-$(RM) $(UTILITIES)/UtilitiesCint.* libUtilities.so
-@echo ' '
-$(RM) $(ACULDATAOBJS) $(ACULDATACPP_DEPS) -$(RM) $(ACULDATAOBJS) $(ACULDATACPP_DEPS)
-$(RM) $(ACULDATA)/AculDataCint.* libAculData.so -$(RM) $(ACULDATA)/AculDataCint.* libAculData.so
-@echo ' ' -@echo ' '
-$(RM) $(ACULCALIBOBJS) $(ACULCALIBCPP_DEPS)
-$(RM) $(ACULCALIB)/AculCalibCint.* libAculCalib.so
-@echo ' '
-$(RM) $(TELOSSOBJS) $(TELOSSCPP_DEPS) -$(RM) $(TELOSSOBJS) $(TELOSSCPP_DEPS)
-$(RM) $(TELOSS)/TELossCint.* libTELoss.so -$(RM) $(TELOSS)/TELossCint.* libTELoss.so
-@echo ' ' -@echo ' '
...@@ -50,17 +66,34 @@ clean: ...@@ -50,17 +66,34 @@ clean:
-@echo ' ' -@echo ' '
# Those *Cint* files below need special treating: # Those *Cint* files below need special treating:
$(UTILITIES)/UtilitiesCint.cpp:
-@echo 'Pre-building UtilitiesCint.cpp and UtilitiesCint.h files'
-rootcint -f $(UTILITIES)/UtilitiesCint.cpp -c -p $(UTILITIES_HEADERS)
-@echo ' '
$(ACULDATA)/AculDataCint.cpp: $(ACULDATA)/AculDataCint.cpp:
-@echo 'Pre-building AculDataCint.cpp and AculDataCint.h files' -@echo 'Pre-building AculDataCint.cpp and AculDataCint.h files'
-rootcint -f $(ACULDATA)/AculDataCint.cpp -c -p $(ACULDATA_HEADERS) -rootcint -f $(ACULDATA)/AculDataCint.cpp -c -p $(ACULDATA_HEADERS)
-@echo ' ' -@echo ' '
$(ACULCALIB)/AculCalibCint.cpp:
-@echo 'Pre-building AculCalibCint.cpp and AculCalibCint.h files'
-rootcint -f $(ACULCALIB)/AculCalibCint.cpp -c -p $(ACULCALIB_HEADERS)
-@echo ' '
$(TELOSS)/TELossCint.cpp: $(TELOSS)/TELossCint.cpp:
-@echo 'Pre-building TELossCint.cpp and TELossCint.h files' -@echo 'Pre-building TELossCint.cpp and TELossCint.h files'
-rootcint -f $(TELOSS)/TELossCint.cpp -c -p $(TELOSS)/TELoss.h $(TELOSS)/linkdef.h -rootcint -f $(TELOSS)/TELossCint.cpp -c -p $(TELOSS)/TELoss.h $(TELOSS)/linkdef.h
-@echo ' ' -@echo ' '
#*.so files #*.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) libAculData.so: libTELoss.so $(ACULDATAOBJS)
@echo 'Building target: $@' @echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker' @echo 'Invoking: GCC C++ Linker'
...@@ -68,6 +101,13 @@ libAculData.so: libTELoss.so $(ACULDATAOBJS) ...@@ -68,6 +101,13 @@ libAculData.so: libTELoss.so $(ACULDATAOBJS)
@echo 'Finished building target: $@' @echo 'Finished building target: $@'
@echo ' ' @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) libTELoss.so: $(TELOSSOBJS)
@echo 'Building target: $@' @echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker' @echo 'Invoking: GCC C++ Linker'
......