Added class for calibration of CsI detectors using beam particles.

AculData/AculCalibCsI.cpp

+#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;
+}

AculData/AculCalibCsI.h

+#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 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();
+};

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

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

macros/calibration_CsI/calibrationSQ13Alpha.cxx

+{
+//	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;
+}

macros/calibration_CsI/calibrationSQ13Trit.cxx

+{
+	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;
+}

macros/calibration_CsI/calibrationSQ23Alpha.cxx

+{
+	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;
+}

macros/calibration_CsI/calibrationSQ23Trit.cxx

+{
+	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");
+
+}

macros/calibration_CsI/parameters/SQ13Alpha.par

+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

macros/calibration_CsI/parameters/SQ13Trit.par

+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

macros/calibration_CsI/parameters/SQ23Alpha.par

+#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

macros/calibration_CsI/parameters/SQ23Trit.par

+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