#include <iostream>
#include <vector>       // std::vector
#include <algorithm>

using namespace std;

void halfVoltPositions680()
{

	TString inFile = "objects/histHalfVolts680.root";

	TString outPeaksFile = "parameters/halfVoltPositions680.par";

	TFile *fr = new TFile(inFile.Data());
	if (fr->IsZombie()) {
		Error("halfVoltPositions680.cxx", "File %s was not open correctly.", inFile.Data());
		return;
	}

	TH1D *hch[64];
	TGraph *gr[64];
	TString hName;

	for (Int_t i = 0; i <64; i++) {
		hName.Form("ch%d", i);
		hch[i] = (TH1D*)fr->Get(hName.Data());
	}



	TString cName;
	TString cTitle;

	TCanvas *cWork[11];
	for (Int_t i = 0; i < 11; i++) {
		cWork[i] = new TCanvas();
		cName.Form("can%d", i);
		cTitle.Form("canvas %d", i);
		cWork[i]->SetName(cName.Data());
		cWork[i]->SetTitle(cTitle.Data());
		cWork[i]->Divide(2, 3);
		cWork[i]->cd(1);
	}

	ofstream calFile;
	calFile.open(outPeaksFile.Data());
	if (calFile.is_open()) {
		Info("calPeaksVoltage.cxx", "Positions of 0.5V peak will be saved in %s.", outPeaksFile.Data());
	} else {
		Error("calPeaksVoltage.cxx", "File %s was not open.", outPeaksFile.Data());
		return;
	}

	Double_t x, y;

	for (Int_t i = 0; i < 64; i++) {
		Int_t cNumber = i/6;
		cWork[cNumber]->cd(i-6*cNumber+1);
		hch[i]->Draw();
		if (i!=42 && i!=21 && i!=11 && i!=41) gr[i] = searchPeaks(hch[i]);
		if (i == 11) {
			gr[i] = searchPeaks(hch[i], 0.2);
		}
		if (i == 21) {
			gr[i] = searchPeaks(hch[i], 0.08);
		}
		if (i == 41) {
			gr[i] = searchPeaks(hch[i], 0.08);
		}
		if (i == 42) {
			gr[i] = searchPeaks(hch[i], 0.08);
		}
		//take a peak corresponding to 0.5 V for 680 pF
		gr[i]->GetPoint(1, x, y);
		calFile << i << "\t" << x << endl;
	}

	calFile.close();

	return;
}

TGraph* searchPeaks(TH1D *hInput, Double_t relHight = 0.05) {

	Double_t voltage[] = {0.5, 1., 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.25, 0.};
//	const Double_t corr1e = 1.18;
//	const Double_t corr1e = 1.;

	hInput->GetXaxis()->SetRangeUser(150., 1000.);

	hInput->ShowPeaks(6., "", relHight);

	TList *functions = hInput->GetListOfFunctions();
	TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");

	const Int_t noPeaks = pm->GetN();
	Double_t *energy = pm->GetX();
	std::vector<Double_t> vEnergy;

	for (Int_t i = 0; i < pm->GetN(); i++) {
		vEnergy.push_back(energy[i]);
	}

	std::sort(vEnergy.begin(), vEnergy.end());

	TGraph *gr = new TGraph(noPeaks+1);
	const Double_t delta_vEnergy = 2*vEnergy[0] - vEnergy[1];

	vEnergy.push_back(delta_vEnergy);
	gr->SetPoint(0, vEnergy[10], voltage[10]/voltage[0]*(vEnergy[0]-delta_vEnergy));
	if (noPeaks == 2) Info("searchPeaks", "%d peaks found in hist %s", noPeaks, hInput->GetName());
	else Error("searchPeaks", "%d peaks found in hist %s", noPeaks, hInput->GetName());
	for (Int_t i = 0; i < noPeaks; i++) {
		if (i>=10) break;
		gr->SetPoint(i+1, vEnergy[i], voltage[i]/voltage[0]*(vEnergy[0]-delta_vEnergy));
	}

	return gr;
}