#include <TROOT.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TGraph.h>

#include <stdint.h>
#include <cstdio>

#include "fits.h"
#include "FOpenDataFile.c"
#include "FOpenCalibFile.c"
#include "FPedestal.c"
#include "FIntFixedPosAllPx.c"

#define n_peaks 2

int poissonanalysis(const char *name, const char *drsname, UInt_t integration_delay, UInt_t integration_size, UInt_t spect_max )
{
//******************************************************************************
//Read a datafile and plot the DRS-calibrated data
//ATTENTION: only works for ROI=1024
// (array indices of the calibration wrong otherwise)
//******************************************************************************

	gROOT->SetStyle("Plain");
//	TFile f("spectra.root","RECREATE");
//-------------------------------------------
//Define the data variables
//-------------------------------------------
	vector<int16_t> data;
	vector<int16_t> data_offset;
	unsigned int data_num;
	size_t data_n;
	UInt_t data_px;
	UInt_t data_roi;
	
//-------------------------------------------
//Get the DRS calibration
//-------------------------------------------
	size_t drs_n;
	vector<float> drs_basemean;
	vector<float> drs_gainmean;
	vector<float> drs_triggeroffsetmean;
	FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
	
//-------------------------------------------
//Check the sizes of the data columns
//-------------------------------------------
//	if(drs_n!=data_n)
//	{
//		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
//		return 1;
//	}
	
//-------------------------------------------
//Open the file
//-------------------------------------------
	fits datafile(name);
	if (!datafile)
	{
	cout << "Couldn't properly open the datafile." << endl;
	return 1;
	}
	FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
	
//-------------------------------------------
//Create the canvas, plot and title
//-------------------------------------------
	char title_spect[500];
	char name_spect[50];
//	std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
	TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
	TH1F* spectrum[data_px];
	for(int i=0; i<data_px; i++) {
		std::sprintf(title_spect,"Spectrum of Px %i, %i-%i (data: %s, DRS: %s)",i,integration_delay,integration_delay+integration_size,name,drsname);
		std::sprintf(name_spect,"spect%i",i);
//		spectrum[i] = new TH1F(name_spect,title_spect,400,-150,650);
		spectrum[i] = new TH1F(name_spect,title_spect,400,-150,spect_max);
		spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
		spectrum[i]->GetYaxis()->SetTitle("Counts");
	}
	
//	TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
//-------------------------------------------
//Define the fit
//-------------------------------------------
//	int fit_start = 20;
//	int fit_end = 200;
//	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp((-x-[7])/[8])",fit_start,fit_end);
//	Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
//	Double_t parread_tmp[5+2*n_peaks];
//	Double_t parread[data_px][5+2*n_peaks];
////	TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
////	Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
//	f_peaks->SetParameters(par_init);
//	f_peaks->SetLineColor(2);
//	
//	TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
//	f_pedestal->SetLineColor(2);
//	f_pedestal->SetLineWidth(1);
//	
//	TF1* f_peak[n_peaks];
//	for(int i=0; i<n_peaks; i++) {
//		f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
//		f_peak[i]->SetLineColor(2);
//		f_peak[i]->SetLineWidth(1);
//	}
	
//-------------------------------------------
//Find the baseline
//-------------------------------------------
//	float pedestal_mean, pedestal_rms;
////	FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
//	//Return values
//	//Value of maximal probability: -2.225 +- 4.3144
//	pedestal_mean = -2.225;
//	pedestal_rms = 4.3144;
	
//-------------------------------------------
//Oscilloscope
//-------------------------------------------
//	float threshold = pedestal_mean+pedestal_rms;
//	float threshold = 2.5;

//***********************
//	UInt_t integration_size = 5;
//	UInt_t integration_delay = 238;

	FIntFixedPosAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, integration_size, integration_delay, spectrum);
	
	canv_spect->cd();
	spectrum[0]->Draw();
	canv_spect->Modified();
	canv_spect->Update();
	
//	Double_t parread_dist[data_px];
//	Double_t parread_x[data_px];
//	for(int i=0; i<data_px; i++) {
//		f_peaks->SetParameters(par_init);
//		if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
////			std::cout << "Fit successful." << std::endl;
//			f_peaks->GetParameters(&parread_tmp[0]);
//			for(int j=0; j<5+2*n_peaks; j++) {
//				parread[i][j] = parread_tmp[j];
//			}
//			parread_dist[i] = parread_tmp[1];
//			parread_x[i] = i;
//		}
//		else {
//			for(int j=0; j<5+2*n_peaks; j++) {
//				parread[i][j] = 0;
//			}
//			parread_dist[i] = 0;
//			parread_x[i] = i;
//		}
//	}
//	
//	canv_gain->cd();
//	TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
//	plot_gain->GetYaxis()->SetRangeUser(0,150);
//	plot_gain->Draw("A*");
//	canv_gain->Modified();
//	canv_gain->Update();
//	
//	canv_spect->cd();
//	canv_spect->SetLogy();
	
//-------------------------------------------
//Draw the data where the gain is out of limits
//-------------------------------------------
//	char temp = 'x';
//	std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
//	for(int i=0; i<data_px; i++) {
//		if((parread_dist[i]<80)||(parread_dist[i]>105)) {
//			std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
//			spectrum[i]->Draw();
//			
//			f_peaks->SetParameters(&parread[i][0]);
//			f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
//			f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
////			f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
////			f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
////			f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
//			f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
//			
//			f_peaks->Draw("same");
//			for(int i=0; i<n_peaks; i++) {
//				f_peak[i]->Draw("same");
//			}
//			f_pedestal->Draw("same");
//			
//			canv_spect->Modified();
//			canv_spect->Update();
//			temp='x';
//			while ((temp!='\n') && (temp!='a'))
//			    cin.get(temp);
//			if(temp=='a') break;
//		}
//	}
//-------------------------------------------
//Draw the data
//-------------------------------------------
	canv_spect->cd();
	Int_t plotspect_nr;
	std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
	std::cin >> plotspect_nr;
	while(plotspect_nr>=0) {
		plotspect_nr = plotspect_nr % data_px;
		spectrum[plotspect_nr]->Draw();
		
//		f_peaks->SetParameters(&parread[plotspect_nr][0]);
//		f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
//		f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
////		f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
////		f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
////		f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
//		f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
//		
//		f_peaks->Draw("same");
//		for(int i=0; i<n_peaks; i++) {
//			f_peak[i]->Draw("same");
//		}
//		f_pedestal->Draw("same");
		
		canv_spect->Modified();
		canv_spect->Update();
		
		std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
		std::cin >> plotspect_nr;
	}
	
//	f->Write();
	return 0;
}