AculCalibScint.cpp

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