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

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

#define HAVE_ZLIB
#include "fits.h"
#include "FOpenDataFile.c"
#include "FOpenCalibFile.c"
#include "FPedestal.c"
#include "FOscilloscopeAllPx.c"

#define n_peaks 2

int gainanalysis(const char *name, const char *drsname, size_t pixelnr)
{
//******************************************************************************
//Read a file with G-APD singles, analyse them and plot the spectra
//This source code is best read with a wide editor (>100 characters per line) or a horizontal scroll bar.
//ATTENTION: this script only works for ROI=1024
// (array indices of the calibration wrong otherwise)
//******************************************************************************

	gROOT->SetStyle("Plain");
	TFile outfile("gainspectra.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 (data: %s, DRS: %s)",i,name,drsname);
		std::sprintf(name_spect,"spect%i",i);
		spectrum[i] = new TH1F(name_spect,title_spect,350,-100,600);
		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};
//	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 (search for peaks and fill integral into the histograms)
//-------------------------------------------
//	float threshold = pedestal_mean+pedestal_rms;
	float threshold = 2.5;

//***********************
	FOscilloscopeAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
	
	canv_spect->cd();
	spectrum[0]->Draw();
	canv_spect->Modified();
	canv_spect->Update();
//***********************
//	fits datafile2("20110804_025.fits");
//	if (!datafile2)
//	{
//	cout << "Couldn't properly open the datafile2." << endl;
//	return 1;
//	}
//	FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px);
//	FOscilloscopeAllPx(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);//	
//	canv_spect->cd();
//	spectrum[0]->Draw();
//	canv_spect->Modified();
//	canv_spect->Update();
////***********************
//	fits datafile3("20110804_026.fits");
//	if (!datafile3)
//	{
//	cout << "Couldn't properly open the datafile3." << endl;
//	return 1;
//	}
//	FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px);
//	FOscilloscopeAllPx(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
//	
//	canv_spect->cd();
//	spectrum[0]->Draw();
//	canv_spect->Modified();
//	canv_spect->Update();
////***********************
//	fits datafile4("20110804_027.fits");
//	if (!datafile4)
//	{
//	cout << "Couldn't properly open the datafile4." << endl;
//	return 1;
//	}
//	FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px);
//	FOscilloscopeAllPx(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
//	
//	canv_spect->cd();
//	spectrum[0]->Draw();
//	canv_spect->Modified();
//	canv_spect->Update();
//***********************
	
//-------------------------------------------
//Fit the histograms and store the parameters into parread
//parread_dist: peak difference (gain) per pixel, parread_x: x-Values for the TGraph
//-------------------------------------------
	Double_t parread[data_px][5+2*n_peaks];
	Double_t parread_tmp[5+2*n_peaks];
	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;
		}
	}
	
//-------------------------------------------
//Plot the gain of all pixel
//-------------------------------------------
	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 spectrum of pixel where the gain is out of limits
//Drawn are the fitted function (f_peaks) and its components (f_peak[i], f_pedestal)
//-------------------------------------------
	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;
		}
	}
	
//-------------------------------------------
//Give the possibility to draw the spectrum of any pixel
//-------------------------------------------
	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;
	}
	outfile.Write();
	outfile.Close();
	return 0;
}