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

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

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

#define n_peaks 4

int peakfinderscope(const char *name, const char *drsname, size_t pixelnr)
{
//******************************************************************************
//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");
	
//-------------------------------------------
//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;
//	}
	
//-------------------------------------------
//Create the canvas, plot and title
//-------------------------------------------
	char title[500];
	std::sprintf(title,"Data: %s, DRS: %s, Px %i",name,drsname,pixelnr);
	TCanvas *canv = new TCanvas( "canv", "Waveform", 200, 10, 700, 500 );
//	TProfile *pulseshape = new TProfile("pulseshape", title, 70, -10.5, 59.5);
	TH2F *pulseshape = new TH2F("pulseshape",title,70,-10.5,59.5,350,-10.5,59.5);
	pulseshape->GetXaxis()->SetTitle("Time after trigger (bins @2 GHz)");
	pulseshape->GetYaxis()->SetTitle("Amplitude (mV)");
	
	char title_base[500];
	std::sprintf(title_base,"All samples of Px %i (data: %s, DRS: %s)",pixelnr,name,drsname);
	
	char title_spect[500];
	std::sprintf(title_spect,"Spectrum of all pixels (data: %s, DRS: %s)",name,drsname);
	TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
	TH1F *spectrum = new TH1F("spectrum", title_spect,700,-200,1200);
	spectrum->GetXaxis()->SetTitle("Pulse size (a.u.)");
	spectrum->GetYaxis()->SetTitle("Counts");
	
//-------------------------------------------
//Define the fit
//-------------------------------------------
	int fit_start = 30;
	int fit_end = 380;
	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-[11])/[12])",fit_start,fit_end);
	Double_t par[5+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,10000,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[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
	f_peaks->SetParameters(par);
	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);
	}
	
//-------------------------------------------
//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);
	
//-------------------------------------------
//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;

//***********************
	FOscilloscope(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
	canv->cd();
	pulseshape->Draw();
	canv->Modified();
	canv->Update();
	
	canv_spect->cd();
	spectrum->Draw();
	canv_spect->Modified();
	canv_spect->Update();
//***********************
	spectrum->Fit(f_peaks,"R+");
	f_peaks->GetParameters(&par[0]);
	
	f_peak[0]->SetParameters(par[2],par[0],par[3]);
	f_peak[1]->SetParameters(par[4],par[0]+par[1],par[5]);
	f_peak[2]->SetParameters(par[6],par[0]+2*par[1],par[7]);
	f_peak[3]->SetParameters(par[8],par[0]+3*par[1],par[9]);
//	f_peak[4]->SetParameters(par[10],par[12],par[11]);
	f_pedestal->SetParameters(par[10],par[11],par[12]);
	
	std::cout << "Crosstalk propability approx. " << ((par[4]*par[5])/(par[2]*par[3]+par[4]*par[5])) << std::endl;
//-------------------------------------------
//Draw the data
//-------------------------------------------
	canv->cd();
	pulseshape->Draw();
	canv->Modified();
	canv->Update();
	
	canv_spect->cd();
	canv_spect->SetLogy();
	spectrum->Draw();
	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();
	
	return 0;
}