From 08178659693a6ca69a2d44a20197135d0a437efe Mon Sep 17 00:00:00 2001 From: Vratislav Chudoba Date: Thu, 20 Oct 2016 12:34:41 +0300 Subject: [PATCH] Added class for calibration of CsI detectors using beam particles. --- AculData/AculCalibCsI.cpp | 655 ++++++++++++++++++ AculData/AculCalibCsI.h | 103 +++ AculData/AculData.mk | 8 +- AculData/linkdef.h | 3 +- .../calibration_CsI/calibrationSQ13Alpha.cxx | 58 ++ .../calibration_CsI/calibrationSQ13Trit.cxx | 37 + .../calibration_CsI/calibrationSQ23Alpha.cxx | 55 ++ .../calibration_CsI/calibrationSQ23Trit.cxx | 46 ++ .../calibration_CsI/parameters/SQ13Alpha.par | 91 +++ .../calibration_CsI/parameters/SQ13Trit.par | 70 ++ .../calibration_CsI/parameters/SQ23Alpha.par | 92 +++ .../calibration_CsI/parameters/SQ23Trit.par | 70 ++ 12 files changed, 1283 insertions(+), 5 deletions(-) create mode 100644 AculData/AculCalibCsI.cpp create mode 100644 AculData/AculCalibCsI.h create mode 100644 macros/calibration_CsI/calibrationSQ13Alpha.cxx create mode 100644 macros/calibration_CsI/calibrationSQ13Trit.cxx create mode 100644 macros/calibration_CsI/calibrationSQ23Alpha.cxx create mode 100644 macros/calibration_CsI/calibrationSQ23Trit.cxx create mode 100644 macros/calibration_CsI/parameters/SQ13Alpha.par create mode 100644 macros/calibration_CsI/parameters/SQ13Trit.par create mode 100644 macros/calibration_CsI/parameters/SQ23Alpha.par create mode 100644 macros/calibration_CsI/parameters/SQ23Trit.par diff --git a/AculData/AculCalibCsI.cpp b/AculData/AculCalibCsI.cpp new file mode 100644 index 0000000..0b3ca75 --- /dev/null +++ b/AculData/AculCalibCsI.cpp @@ -0,0 +1,655 @@ +#include "AculCalibCsI.h" + +ClassImp(AculCalibCsI); + +AculCalibCsI::AculCalibCsI() : fr("TFile", 5), tr("TTree", 5), cutsCol("TCutG"), fA(16), fB(16) { + + printf("AculCalibCsI::Default constructor called.\n"); + + +} + +AculCalibCsI::AculCalibCsI(const char* parfile) : fr("TFile", 5), tr("TTree", 5), cutsCol("TCutG", 10), fA(16), fB(16) { + printf("AculCalibCsI::Constructor called.\n"); + + SetParFile(parfile); + SetPars(); + cout << "nofiles: " << nofiles << endl; + OpenTrees(); + LoadCuts(); + +// fr.Print(); +// fr.At(0); + + +} + +AculCalibCsI::~AculCalibCsI() { + + printf("AculCalibCsI::Destructor called.\n"); + +} + +void AculCalibCsI::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 AculCalibCsI::LoadCuts() { + + if (cutsFileName == "") { + printf("AculCalibCsI::LoadCuts: Name of file (*.root) with cuts was not provided.\n"); + printf("AculCalibCsI::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 AculCalibCsI::SetParFile(const char* parfile) { + + fParFileName = parfile; + return; +} + +void AculCalibCsI::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 AculCalibCsI::PrintPeakRanges() { + for (Int_t i = 0; i < 16; i++) { + printf("%d\t%d\t%d\n", i, peakMin[0][i], peakMax[0][i]); + } +} + +void AculCalibCsI::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 AculCalibCsI::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 AculCalibCsI::PrintCuts() { + +// printf("AculCalibCsI::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("AculCalibCsI::PrintCuts: End of function.\n"); +*/ + return; +} + +Double_t AculCalibCsI::GetA(Int_t i) { + if (i >= fA.GetSize()) //if i >= number of array element + { + return 0.; + } + return fA[i]; +} + +Double_t AculCalibCsI::GetB(Int_t i) { + if (i >= fB.GetSize()) //if i >= number of array element + { + return 0.; + } + return fB[i]; +} + +void AculCalibCsI::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("AculCalibCsI::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 AculCalibCsI::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("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 < energyPoints; j++) { + if ( cutsCol.At(j) ) { + curCutG = (TCutG*)cutsCol.At(j); + curCutG->Draw("same"); +// printf("AculCalibCsI::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("AculCalibCsI::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("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j); +// } + } + canvas->Update(); + } +} + +void AculCalibCsI::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("AculCalibCsI::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 AculCalibCsI::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("AculCalibCsI::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 AculCalibCsI::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 AculCalibCsI::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 AculCalibCsI::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" << 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 AculCalibCsI::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("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 AculCalibCsI::ReadClbParameters(const char* filename) { + + std::ifstream infile(filename); + if ( !infile.is_open() ) { + printf("AculCalibCsI::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 AculCalibCsI::DrawClbGraphs(const char* filename, const char* graphname, TCanvas *canvas) { + + printf("AculCalibCsI::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 AculCalibCsI::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("AculCalibCsI::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(); + } + + +} + +void AculCalibCsI::SetPars() { + + std::ifstream infile(fParFileName.Data()); + if ( !infile.is_open() ) { + printf("AculCalibCsI::ReadParFile: File %s was not open.\n", fParFileName.Data()); + 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("AculCalibCsI::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; +} diff --git a/AculData/AculCalibCsI.h b/AculData/AculCalibCsI.h new file mode 100644 index 0000000..c3a0add --- /dev/null +++ b/AculData/AculCalibCsI.h @@ -0,0 +1,103 @@ +#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 AculCalibCsI { + +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; + + TString fParFileName; + +public: +// AculCalibCsI() : a(0), b(0), c(0), p(0){}; + AculCalibCsI(); + AculCalibCsI(const char* parfile); + virtual ~AculCalibCsI(); + // Define the class for the cint dictionary + ClassDef (AculCalibCsI,1); + + + void OpenTrees(); + void LoadCuts(); + + void SetParFile(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); + Double_t GetB(Int_t i); + +private: + + void SetPars(); +}; diff --git a/AculData/AculData.mk b/AculData/AculData.mk index ad49fcb..0d57c29 100755 --- a/AculData/AculData.mk +++ b/AculData/AculData.mk @@ -7,20 +7,20 @@ ACULDATALIBS := -lCore -lCint -lRIO -lTree -lNet -lThread -lHist -lMatrix -lMath # Add inputs and outputs from these tool invocations to the build variables ACULDATA_HEADERS += \ $(ACULDATA)/AculCalibration.h \ -$(ACULDATA)/ConfigDictionary.h \ +$(ACULDATA)/AculCalibCsI.h \ $(ACULDATA)/linkdef.h ACULDATACPP_SRCS += \ $(ACULDATA)/AculCalibration.cpp \ -$(ACULDATA)/ConfigDictionary.cpp \ +$(ACULDATA)/AculCalibCsI.cpp \ $(ACULDATA)/AculDataCint.cpp ACULDATAOBJS += \ $(ACULDATA)/AculCalibration.o \ -$(ACULDATA)/ConfigDictionary.o \ +$(ACULDATA)/AculCalibCsI.o \ $(ACULDATA)/AculDataCint.o ACULDATACPP_DEPS += \ $(ACULDATA)/AculCalibration.d \ -$(ACULDATA)/ConfigDictionary.d \ +$(ACULDATA)/AculCalibCsI.d \ $(ACULDATA)/AculDataCint.d \ No newline at end of file diff --git a/AculData/linkdef.h b/AculData/linkdef.h index 24524d2..42dc2a3 100755 --- a/AculData/linkdef.h +++ b/AculData/linkdef.h @@ -6,7 +6,8 @@ //#pragma link C++ class AculRaw; //#pragma link C++ class AculConvert; #pragma link C++ class AculCalibration; -#pragma link C++ class ConfigDictionary; +#pragma link C++ class AculCalibCsI; +//#pragma link C++ class ConfigDictionary; #endif diff --git a/macros/calibration_CsI/calibrationSQ13Alpha.cxx b/macros/calibration_CsI/calibrationSQ13Alpha.cxx new file mode 100644 index 0000000..3a337aa --- /dev/null +++ b/macros/calibration_CsI/calibrationSQ13Alpha.cxx @@ -0,0 +1,58 @@ +{ +// gSystem->Load("CsICalib_C.so"); + gSystem->Load("../../libAculData.so"); + + //CsICalib cal; + AculCalibCsI cal("parameters/SQ13Alpha.par"); +// CsICalib cal("SQ13AlphaAlt.par"); + +// cal.OpenSQ13AlphaTrees(); +// cal.LoadCutsSQ13A(); + + TCanvas *c1 = new TCanvas("c1", "Plain"); + TCanvas *c2 = new TCanvas("c2", "Plain"); + TCanvas *c3 = new TCanvas("c3", "Plain"); + TCanvas *c4 = new TCanvas("c4", "Plain"); + + cal.PrintTrees(); + cal.PrintFiles(); + cal.PrintCuts(); + +// cal.ReadParFile("ranges.dat"); +// cal.PrintParameters(); + + + +// cal.DrawVariable("SQ13", 1, c1); +//return; +// cal.DrawVariable("TDC", 2, c2); +// cal.DrawVariable("QDC", 3, c3); + +// cal.DrawBeam(c4, 4, "SQ13"); +//return; +// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16"); +// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp"); +//return; +// cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp"); +// cal.DrawVariableCut("SQ13", 2, c3, "cutSQ13Alpha26", "cutSQ13Alpha26Amp"); +// cal.DrawVariableCut("SQ13", 3, c4, "cutSQ13Alpha30", "cutSQ13Alpha30Amp"); +//return; + + cal.GetPeakMean("SQ13", 0, 0, c1, "cutSQ13Alpha16Amp", 256); + cal.GetPeakMean("SQ13", 1, 1, c2, "cutSQ13Alpha21Amp", 256); + cal.GetPeakMean("SQ13", 2, 2, c3, "cutSQ13Alpha26Amp", 256); + cal.GetPeakMean("SQ13", 3, 3, c4, "cutSQ13Alpha30Amp", 256); + + TCanvas *cCal = new TCanvas("cCal", "calibration Alpha"); + cal.Calibrate(cCal); + +// cal.SaveClbGraphs("gSQ13Alpha.root", "RECREATE"); +// cal.SaveClbGraphs("gSQ13Alpha2points.root", "RECREATE"); +// cal.SaveClbGraphs("gSQ13AlphaAlt.root", "RECREATE"); + +// cal.WriteClbParameters("calSQ13Alpha.clb"); +// cal.WriteClbParameters("calSQ13Alpha2points.clb"); +// cal.WriteClbParameters("calSQ13AlphaAlt.clb"); + +return; +} diff --git a/macros/calibration_CsI/calibrationSQ13Trit.cxx b/macros/calibration_CsI/calibrationSQ13Trit.cxx new file mode 100644 index 0000000..e7bd12f --- /dev/null +++ b/macros/calibration_CsI/calibrationSQ13Trit.cxx @@ -0,0 +1,37 @@ +{ + gSystem->Load("CsICalib_C.so"); + + //CsICalib cal; + CsICalib cal("SQ13Trit.par"); + +// cal.OpenSQ13TritTrees(); +// cal.LoadCutsSQ13T(); + + TCanvas *c1 = new TCanvas("c1", "Plain"); + TCanvas *c2 = new TCanvas("c2", "Plain"); + TCanvas *c3 = new TCanvas("c3", "Plain"); + +// cal.PrintTrees(); +// cal.PrintFiles(); +// cal.PrintCuts(); + + cal.PrintParameters(); + +// cal.DrawBeam(c3, 3); + +// cal.DrawVariableCut("SQ13", 0, c1, "cutSQ13Trit"); +// cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Trit19", "cutSQ13Trit19Amp"); +// cal.DrawVariableCut("SQ13", 2, c3, "cutSQ13Trit24", "cutSQ13Trit24Amp"); + + cal.GetPeakMean("SQ13", 0, 0, c1, "cutSQ13Trit14Amp", 256); + cal.GetPeakMean("SQ13", 1, 1, c2, "cutSQ13Trit19Amp", 256); + cal.GetPeakMean("SQ13", 2, 2, c3, "cutSQ13Trit24Amp", 256); + + + TCanvas *cCal = new TCanvas("cCal", "calibration SQ13 Triton"); + cal.Calibrate(cCal); + cal.SaveClbGraphs("gSQ13Trit.root", "RECREATE"); + cal.WriteClbParameters("calSQ13Trit.clb"); + +return; +} diff --git a/macros/calibration_CsI/calibrationSQ23Alpha.cxx b/macros/calibration_CsI/calibrationSQ23Alpha.cxx new file mode 100644 index 0000000..06fb76c --- /dev/null +++ b/macros/calibration_CsI/calibrationSQ23Alpha.cxx @@ -0,0 +1,55 @@ +{ + gSystem->Load("CsICalib_C.so"); + +// CsICalib cal; + CsICalib cal("parameters/SQ23Alpha.par"); + + +// cal.PrintCuts(); + +//return; + + TCanvas *c1 = new TCanvas("c1", "Plain"); + TCanvas *c2 = new TCanvas("c2", "Plain"); + TCanvas *c3 = new TCanvas("c3", "Plain"); + TCanvas *c4 = new TCanvas("c4", "Plain"); + + cal.DrawBeam(c4, 4, "SQ23"); + +return; + +// cal.ReadParFile("SQ23Alpha.par"); +// cal.OpenTrees(); + +// cal.ReadClbParameters("calSQ23Alpha.clb"); +// cal.DrawEnergyDeposite("SQ23", c4, 3, "same"); + +//return; + +// cal.DrawVariable("SQ23", 1, c1); +// cal.DrawVariable("TDC", 2, c2, 0, 700); +// cal.DrawVariable("QDC", 3, c3); +//return; + +// cal.DrawVariableCut("SQ23", 0, c1, "cutAlpha16", "cutAlpha16Amp"); +// cal.DrawVariableCut("SQ23", 0, c2, "cutAlpha20", "cutAlpha20Amp"); +// cal.DrawVariableCut("SQ23", 0, c3, "cutAlpha25", "cutAlpha25Amp"); +// cal.DrawVariableCut("SQ23", 0, c4, "cutAlpha30", "cutAlpha30Amp"); + +// cal.GetPeakMean("SQ23", 0, c1, "cutAlpha16Amp", minAlpha16, maxAlpha16); +// cal.GetPeakMean("SQ23", 1, c2, "cutAlpha20Amp", minAlpha20, maxAlpha20); +// cal.GetPeakMean("SQ23", 2, c3, "cutAlpha25Amp", minAlpha25, maxAlpha25); +// cal.GetPeakMean("SQ23", 3, c4, "cutAlpha30Amp", minAlpha30, maxAlpha30); + + cal.GetPeakMean("SQ23", 0, 0, c1, "cutAlpha16Amp"); + cal.GetPeakMean("SQ23", 1, 1, c2, "cutAlpha20Amp"); + cal.GetPeakMean("SQ23", 2, 2, c3, "cutAlpha25Amp"); + cal.GetPeakMean("SQ23", 3, 3, c4, "cutAlpha30Amp"); + + TCanvas *cCal = new TCanvas("cCal", "calibration Alpha"); + cal.Calibrate(cCal); +// cal.SaveClbGraphs("gSQ23Alpha.root", "RECREATE"); +// cal.WriteClbParameters("calSQ23Alpha.clb"); + +return; +} diff --git a/macros/calibration_CsI/calibrationSQ23Trit.cxx b/macros/calibration_CsI/calibrationSQ23Trit.cxx new file mode 100644 index 0000000..01045fe --- /dev/null +++ b/macros/calibration_CsI/calibrationSQ23Trit.cxx @@ -0,0 +1,46 @@ +{ + gSystem->Load("CsICalib_C.so"); + + CsICalib cal("SQ23Trit.par"); + + TCanvas *c1 = new TCanvas("c1", "Plain"); + TCanvas *c2 = new TCanvas("c2", "Plain"); + TCanvas *c3 = new TCanvas("c3", "Plain"); + +// cal.LoadCutsSQ23T(); + cal.PrintCuts(); + +// cal.PrintTrees(); +// cal.PrintFiles(); +// cal.PrintCuts(); +// cal.ReadPeakRanges(); +// cal.PrintPeakRanges(); +// cal.ReadParFile("ranges.dat"); +// cal.PrintParameters(); +// return; + + +// cal.DrawVariable("SQ23", 1, c1); +// cal.DrawVariable("TDC", 2, c2); +// cal.DrawVariable("QDC", 3, c3); + +// cal.DrawVariableCut("SQ23", 0, c1, "cutAlpha16", "cutAlpha16Amp"); +// cal.DrawVariableCut("SQ23", 0, c2, "cutAlpha20", "cutAlpha20Amp"); +// cal.DrawVariableCut("SQ23", 0, c3, "cutAlpha25", "cutAlpha25Amp"); +// cal.DrawVariableCut("SQ23", 0, c4, "cutAlpha30", "cutAlpha30Amp"); + +// cal.GetPeakMean("SQ23", 0, c1, "cutAlpha16Amp", minAlpha16, maxAlpha16); +// cal.GetPeakMean("SQ23", 1, c2, "cutAlpha20Amp", minAlpha20, maxAlpha20); +// cal.GetPeakMean("SQ23", 2, c3, "cutAlpha25Amp", minAlpha25, maxAlpha25); +// cal.GetPeakMean("SQ23", 3, c4, "cutAlpha30Amp", minAlpha30, maxAlpha30); + + cal.GetPeakMean("SQ23", 0, 0, c1, "cutTrit13Amp"); + cal.GetPeakMean("SQ23", 1, 1, c2, "cutTrit20Amp"); + cal.GetPeakMean("SQ23", 2, 2, c3, "cutTrit24Amp"); + + TCanvas *cCal = new TCanvas("cCal", "calibration SQ23 Alpha"); + cal.Calibrate(cCal); + cal.SaveClbGraphs("gSQ23Trit.root", "RECREATE"); + cal.WriteClbParameters("calSQ23Trit.clb"); + +} diff --git a/macros/calibration_CsI/parameters/SQ13Alpha.par b/macros/calibration_CsI/parameters/SQ13Alpha.par new file mode 100644 index 0000000..2a1dbbd --- /dev/null +++ b/macros/calibration_CsI/parameters/SQ13Alpha.par @@ -0,0 +1,91 @@ +energies 4 #number of energy points +files 4 +data/csi_13_Ealpha16.root +data/csi_13_Ealpha21.root +data/csi_13_Ealpha26.root +data/csi_13_Ealpha30.root +#keyword; filename; number of cuts in file +cutFile cuts/cutsSQ13Alpha.root 8 +cutsSQ13Alpha16 +cutSQ13Alpha21 +cutSQ13Alpha26 +cutSQ13Alpha30 +cutSQ13Alpha16Amp +cutSQ13Alpha21Amp +cutSQ13Alpha26Amp +cutSQ13Alpha30Amp +#detector particle +detector SQ13 Alpha +#energy [MeV/A] +energy 66.542 +#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 +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 +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 +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 diff --git a/macros/calibration_CsI/parameters/SQ13Trit.par b/macros/calibration_CsI/parameters/SQ13Trit.par new file mode 100644 index 0000000..894451c --- /dev/null +++ b/macros/calibration_CsI/parameters/SQ13Trit.par @@ -0,0 +1,70 @@ +energies 3 #number of energy points +files 3 +csi_13_Etrit14.root +csi_13_Etrit19.root +csi_13_Etrit24.root +#keyword; filename; number of cuts in file +cutFile cutsSQ13Trit.root 6 +cutSQ13Trit +cutSQ13Trit19 +cutSQ13Trit24 +cutSQ13Trit14Amp +cutSQ13Trit19Amp +cutSQ13Trit24Amp +#detector particle +detector SQ13 Triton +#energy [MeV/A] +energy 43.636 #14 AMeV +#channel minimum maximum +0 300 500 +1 400 550 +2 400 500 +3 500 700 +4 450 550 +5 350 450 +6 300 400 +7 500 600 +8 400 550 +9 300 400 +10 450 600 +11 400 600 +12 350 450 +13 550 700 +14 200 300 +15 400 500 +#CsI2, energy ~ 19 MeV/A, Triton +energy 61.52 +0 200 500 +1 650 750 +2 580 680 +3 850 1000 +4 650 800 +5 500 600 +6 450 550 +7 750 850 +8 600 750 +9 450 550 +10 650 850 +11 650 750 +12 500 600 +13 800 1000 +14 250 350 +15 590 700 +#CsI2, energy ~ 24 MeV/A, Triton +energy 69.382 +0 600 700 +1 800 950 +2 700 850 +3 400 600 +4 800 1000 +5 600 700 +6 550 650 +7 950 1050 +8 780 900 +9 550 650 +10 850 1000 +11 800 900 +12 650 740 +13 1050 1200 +14 350 450 +15 740 840 diff --git a/macros/calibration_CsI/parameters/SQ23Alpha.par b/macros/calibration_CsI/parameters/SQ23Alpha.par new file mode 100644 index 0000000..97784cc --- /dev/null +++ b/macros/calibration_CsI/parameters/SQ23Alpha.par @@ -0,0 +1,92 @@ +#keywords: detector, energies, files, cutFile +energies 4 #number of energy points +files 4 +data/csi2_beam_Ealpha16all.root +data/csi2_beam_Ealpha20all.root +data/csi2_beam_Ealpha25all.root +data/csi2_beam_Ealpha30all.root +#keyword; filename; number of cuts in file +cutFile cuts/beamCuts.root 8 +cutAlpha16 +cutAlpha20 +cutAlpha25 +cutAlpha30 +cutAlpha16Amp +cutAlpha20Amp +cutAlpha25Amp +cutAlpha30Amp +#detector particle +detector SQ23 Alpha +#energy [MeV/A] +energy 66.508 +#channel minimum maximum +0 650 830 +1 1060 1200 +2 1040 1180 +3 1080 1200 +4 1060 1200 +5 1250 1460 +6 1100 1280 +7 1200 1400 +8 1100 1250 +9 1050 1200 +10 1050 1220 +11 1100 1280 +12 660 750 +13 1150 1300 +14 1400 1600 +15 1250 1420 +#CsI2, energy ~ 20 MeV/A, Alpha +energy 85.236 +0 900 1050 +1 1400 1550 +2 1350 1500 +3 1400 1550 +4 1400 1550 +5 1650 1900 +6 1450 1650 +7 1600 1800 +8 1400 1650 +9 1350 1550 +10 1400 1600 +11 1400 1650 +12 800 1100 +13 1450 1650 +14 1800 2050 +15 1650 1850 +#CsI2, energy ~ 25 MeV/A, Alpha +energy 104.532 +0 1100 1250 +1 1730 1900 +2 1680 1850 +3 1600 2000 +4 1720 1900 +5 2100 2300 +6 1800 2000 +7 2000 2200 +8 1800 2000 +9 1700 1900 +10 1750 1900 +11 1800 2000 +12 1050 1200 +13 1850 2050 +14 2300 2500 +15 2030 2300 +#CsI2, energy ~ 30 MeV/A, Alpha +energy 119.42 +0 1200 1400 +1 2000 2150 +2 1900 2050 +3 1900 2200 +4 1950 2150 +5 2400 2600 +6 2050 2250 +7 2280 2500 +8 2060 2200 +9 1900 2100 +10 1980 2150 +11 2080 2250 +12 1100 1400 +13 2100 2300 +14 2600 2850 +15 2380 2550 diff --git a/macros/calibration_CsI/parameters/SQ23Trit.par b/macros/calibration_CsI/parameters/SQ23Trit.par new file mode 100644 index 0000000..5fe9743 --- /dev/null +++ b/macros/calibration_CsI/parameters/SQ23Trit.par @@ -0,0 +1,70 @@ +energies 3 #number of energy points +files 3 +csi2_beam_Ealpha30all.root +csi2_Etrit20all.root +csi2_Etrit24all.root +#keyword; filename; number of cuts in file +cutFile beamCuts_CsI2_trit.root 6 +cutTrit13 +cutTrit20 +cutTrit24 +cutTrit13Amp +cutTrit20Amp +cutTrit24Amp +#detector particle +detector SQ23 Tritium +#CsI2, energy ~ 13 MeV/A, Tritium +energy 39.102 +#channel minimum maximum +0 450 520 +1 670 750 +2 650 750 +3 650 800 +4 700 750 +5 800 900 +6 700 800 +7 780 850 +8 700 800 +9 650 750 +10 680 750 +11 700 800 +12 440 500 +13 740 800 +14 880 990 +15 800 900 +#CsI2, energy ~ 20 MeV/A, Tritium +energy 61.17 +0 700 750 +1 1040 1100 +2 1000 1100 +3 1000 1150 +4 1050 1150 +5 1250 1350 +6 1100 1200 +7 1200 1300 +8 1050 1200 +9 1000 1100 +10 1040 1150 +11 1050 1200 +12 650 700 +13 1100 1200 +14 1350 1450 +15 1000 1400 +#CsI2, energy ~ 24 MeV/A, Tritium +energy 73.272 +0 800 900 +1 1200 1300 +2 1200 1300 +3 700 1300 +4 1200 1350 +5 1450 1800 +6 1250 1400 +7 1400 1500 +8 1250 1400 +9 1200 1300 +10 1200 1350 +11 1300 1450 +12 700 850 +13 1300 1400 +14 1600 1800 +15 1450 1550 -- 2.18.1