#include "TH1.h"
#include "TF1.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>

//*******************************
#include <TGClient.h>
#include <TCanvas.h>
#include <TRandom.h>
#include <TGButton.h>
#include <TGFrame.h>
#include <TRootEmbeddedCanvas.h>
#include <TMath.h>

//********************************
using namespace std;
using std::ofstream;

const Int_t max_peak_num = 10;

//**********************************
class MyMainFrame : public TGMainFrame {
private:
   TRootEmbeddedCanvas *fEcanvas;
   TRootEmbeddedCanvas *fEcanvas2;
   fEcanvas2->ToggleEventStatus();
   fEcanvas2->ToggleToolBar();
   TGTextEntry *bt_spe_name;
   TGCheckButton *bt_logscale;
   TGNumberEntry *bt_left;
   TGNumberEntry *bt_right;
   TGNumberEntry *bt_left2;
   TGNumberEntry *bt_right2;
   TGCheckButton *bt_logscale2;
   TGCheckButton *bt_offset;
   TGCheckButton *bt_offset2;
   TGTextEntry *bt_peak_pos;
   TGTextEntry *bt_sig_par;
   TGNumberEntry *bt_Esum; 
   TGNumberEntry *bt_bkg_deg; 
   TH1F *hspe;
   TH1F *hspe1;
   TH1F *hfit;
public:
   MyMainFrame(const TGWindow *p, UInt_t w, UInt_t h);
   virtual ~MyMainFrame();
   void DoDraw();
   void DoDraw2();
};

MyMainFrame::MyMainFrame(const TGWindow *p, UInt_t w, UInt_t h) : TGMainFrame(p, w, h)
{
   // Great a 1st horizontal frame
   TGHorizontalFrame *h1frame = new TGHorizontalFrame(this, 1000, 400);
   TGVerticalFrame *v11frame = new TGVerticalFrame(h1frame, 200, 400, kFixedWidth | kFixedHeight);
   TGVerticalFrame *v12frame = new TGVerticalFrame(h1frame, 600, 400, kFixedWidth | kFixedHeight);
   TGVerticalFrame *v13frame = new TGVerticalFrame(h1frame, 200, 400, kFixedWidth | kFixedHeight);
   
   bt_spe_name = new TGTextEntry(v11frame,"GE5860.dat");
   TGTextButton *bt_load = new TGTextButton(v11frame, "Load spectrum");
   TGTextButton *bt_exit = new TGTextButton(v11frame, "Exit", "gApplication->Terminate(0)");
   bt_logscale = new TGCheckButton(v11frame, "Logscale",0);
   bt_offset = new TGCheckButton(v11frame, "Offset",0);
   bt_left = new TGNumberEntry(v11frame, 0., 0., 16000.);
   bt_right = new TGNumberEntry(v11frame, 8192., 0., 16000.);
   TGLabel *label_sig = new TGLabel(v11frame, "Sigma parameters");
   bt_sig_par = new TGTextEntry(v11frame,"0.000114 1.9 0.00028 1.5"); 
   TGLabel *label_sum = new TGLabel(v11frame, "Energy sum (channel)");
   bt_Esum = new TGNumberEntry(v11frame,5830.);
   
   fEcanvas = new TRootEmbeddedCanvas("ECanvas", v12frame, 200, 200);
   TGLabel *label1 = new TGLabel(v11frame, "Left:");
   TGLabel *label2 = new TGLabel(v11frame, "Right:");

   v11frame->AddFrame(bt_spe_name, new TGLayoutHints(kLHintsExpandX, 4, 4, 0, 0));
   v11frame->AddFrame(bt_load, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_exit, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_logscale, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_offset, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(label1, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_left, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(label2, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_right, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(label_sig, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_sig_par, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(label_sum, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   v11frame->AddFrame(bt_Esum, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
   
   v12frame->AddFrame(fEcanvas, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));

   h1frame->AddFrame(v11frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
   h1frame->AddFrame(v12frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
   h1frame->AddFrame(v13frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
   
   AddFrame(h1frame, new TGLayoutHints(kLHintsExpandX| kLHintsExpandY, 10, 10, 10, 10));

   
   bt_load->Connect("Clicked()", "MyMainFrame", this, "DoDraw()");
   bt_logscale->Connect("Clicked()", "MyMainFrame", this, "DoDraw()");
   bt_offset->Connect("Clicked()", "MyMainFrame", this, "DoDraw()");
	bt_left->Connect("ValueSet(Long_t)", "MyMainFrame", this, "DoDraw()");
	bt_left->GetNumberEntry()->Connect("ReturnPressed()","MyMainFrame", this, "DoDraw()");
	bt_right->Connect("ValueSet(Long_t)", "MyMainFrame", this, "DoDraw()");
	bt_right->GetNumberEntry()->Connect("ReturnPressed()","MyMainFrame", this, "DoDraw()");
	

	// Second parts
	TGHorizontalFrame *h2frame = new TGHorizontalFrame(this, 1000, 400);
	TGVerticalFrame *v21frame = new TGVerticalFrame(h2frame, 200, 400, kFixedWidth | kFixedHeight);
    TGVerticalFrame *v22frame = new TGVerticalFrame(h2frame, 600, 400, kFixedWidth | kFixedHeight);
	TGVerticalFrame *v23frame = new TGVerticalFrame(h2frame, 200, 400, kFixedWidth | kFixedHeight);
	
	TGLabel *label3 = new TGLabel(v21frame, "Fitting Section");
	v21frame->AddFrame(label3, new TGLayoutHints(kLHintsExpandX, 4, 4, 0, 0));
	
	bt_logscale2 = new TGCheckButton(v21frame, "Logscale",0);
	bt_offset2 = new TGCheckButton(v21frame, "Offset",0);
	bt_left2 = new TGNumberEntry(v21frame, 570., 0., 16000.);
	bt_right2 = new TGNumberEntry(v21frame, 640., 0., 16000.);
	TGTextButton *bt_fit = new TGTextButton(v21frame, "Fitting");
	TGLabel *label5 = new TGLabel(v21frame, "Peak positions:"); //Devided by space
	bt_peak_pos = new TGTextEntry(v21frame, "603.");
	TGLabel *label_bkg = new TGLabel(v21frame, "Bkg poly degree (<=2)"); 
	bt_bkg_deg = new TGNumberEntry(v21frame,1,0,2);
	
		
	v21frame->AddFrame(bt_logscale2, new TGLayoutHints(kLHintsExpandX, 4, 4, 0, 0));
	v21frame->AddFrame(bt_offset2, new TGLayoutHints(kLHintsExpandX, 4, 4, 4, 0));
	v21frame->AddFrame(bt_left2, new TGLayoutHints(kLHintsExpandX, 4, 4, 4, 0));
	v21frame->AddFrame(bt_right2, new TGLayoutHints(kLHintsExpandX, 4, 4, 4, 0));
	v21frame->AddFrame(label5, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
	v21frame->AddFrame(bt_peak_pos, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));

	v21frame->AddFrame(label_bkg, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
	v21frame->AddFrame(bt_bkg_deg, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
	v21frame->AddFrame(bt_fit, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
	
	
	fEcanvas2 = new TRootEmbeddedCanvas("ECanvas", v22frame, 200, 200);
	v22frame->AddFrame(fEcanvas2, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
	
	h2frame->AddFrame(v21frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
    h2frame->AddFrame(v22frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
	h2frame->AddFrame(v23frame, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 0, 0, 0, 0));
	
	
   AddFrame(h2frame, new TGLayoutHints(kLHintsExpandX| kLHintsExpandY, 10, 10, 10, 10));

   bt_logscale2->Connect("Clicked()", "MyMainFrame", this, "DoDraw2()");
   bt_offset2->Connect("Clicked()", "MyMainFrame", this, "DoDraw2()");
   bt_left2->Connect("ValueSet(Long_t)", "MyMainFrame", this, "DoDraw2()");
	bt_left2->GetNumberEntry()->Connect("ReturnPressed()","MyMainFrame", this, "DoDraw2()");
	bt_right2->Connect("ValueSet(Long_t)", "MyMainFrame", this, "DoDraw2()");
	bt_right2->GetNumberEntry()->Connect("ReturnPressed()","MyMainFrame", this, "DoDraw2()");
	bt_right2->GetNumberEntry()->Connect("ReturnPressed()","MyMainFrame", this, "DoDraw2()");
	bt_fit->Connect("Clicked()", "MyMainFrame", this, "DoFit()");
	
   // Set a name to the main frame
   SetWindowName("Fitting peaks");
   // Map all subwindows of main frame
   MapSubwindows();
   // Initialize the layout algorithm
   Resize(GetDefaultSize());
   // Map main frame
   MapWindow();
   
   
}

void MyMainFrame::DoDraw()
{
	hspe = loadspe(bt_spe_name->GetText());
	hspe1 = hspe->Clone();
   TCanvas *fCanvas = fEcanvas->GetCanvas();
   fCanvas->cd();
   // Check offset
	if (bt_offset->IsOn())
	{
		if (hspe->GetMinimum()<0) 
		{
		Double_t addval = hspe->GetMinimum();
		//cout<<"Value = "<<addval<<endl;
		for (int i = 1; i<=16084; i++) hspe1->SetBinContent(i, hspe->GetBinContent(i) - addval + 1);
		}  	
	}
   
// Check logscale status
	if (bt_logscale->IsOn() == true) fCanvas->SetLogy();	
	if (bt_logscale->IsOn() == false) fCanvas->SetLogy(0); 

// Check range
	hspe1->GetXaxis()->SetRange(bt_left->GetNumber(),bt_right->GetNumber()); 
   hspe1->Draw();
   fCanvas->Update();
}

void MyMainFrame::DoDraw2()
{
	MyMainFrame::DoDraw();
	TCanvas *fCanvas = fEcanvas2->GetCanvas();
	fCanvas->cd();
	
	
	hfit = hspe->Clone();
	
	// Check logscale status
	if (bt_logscale2->IsOn() == true) fCanvas->SetLogy();
	if (bt_logscale2->IsOn() == false) fCanvas->SetLogy(0); 
	
	// Check offset
	if (bt_offset2->IsOn())
	{
		if (hspe->GetMinimum()<0) 
		{
		Double_t addval = hspe->GetMinimum();
		//cout<<"Value = "<<addval<<endl;
		for (int i = 1; i<=16084; i++) hfit->SetBinContent(i, hspe->GetBinContent(i) - addval + 1);
		}  	
	}
	
	
	hfit->GetXaxis()->SetRange(bt_left2->GetNumber(), bt_right2->GetNumber());
	hfit->SetTitle("Fitting");
	hfit->Draw();
	
	fCanvas->Update();
}

Double_t fit_function(Double_t *x, Double_t *par) //max_peak_num
{
	Double_t ans = 0.; 
	ans = ans + par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
	for (int i=0; i<max_peak_num; i++) 
	{
		ans = ans + par[3+3*i]*TMath::Exp(-2.0*(x[0]-par[4+3*i])*(x[0]-par[4+3*i])/par[5+3*i]/par[5+3*i])/par[5+3*i]/sqrt(3.14/2);
	}
	return ans;
}

Double_t fit_function2(Double_t *x, Double_t *par) //20*4+3 = 83 parameters; two gauss
{
	Double_t ans = 0; 
	ans = ans + par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
	for (int i=0; i<max_peak_num; i++) 
	{
		ans = ans + 0.5*par[3+4*i]*(TMath::Exp(-2.0*(x[0]-par[4+4*i])*(x[0]-par[4+4*i])/par[5+4*i]/par[5+4*i])/par[5+4*i]/sqrt(3.14/2) + TMath::Exp(-2.0*(x[0]-par[4+4*i])*(x[0]-par[4+4*i])/par[6+4*i]/par[6+4*i])/par[6+4*i]/sqrt(3.14/2));
	}
	return ans;
}

Double_t cal_sig(double sp[4], double E1, double E2)
{
	double tmp1 = sp[0]*E1+sp[1]; 
	double tmp2 = sp[2]*E2+sp[3];
	double tmp3 = sp[0]*E2+sp[1];
	double tmp4 = sp[2]*E1+sp[3];
	return 0.5*(sqrt(tmp1*tmp1*tmp2*tmp2/(tmp1*tmp1+tmp2*tmp2))+sqrt(tmp3*tmp3*tmp4*tmp4/(tmp3*tmp3+tmp4*tmp4)));
}

Double_t cal_sig1(double sp[4], double E1, double E2)
{
	double tmp1 = sp[0]*E1+sp[1]; 
	double tmp2 = sp[2]*E2+sp[3];
	return sqrt(tmp1*tmp1*tmp2*tmp2/(tmp1*tmp1+tmp2*tmp2));
}

Double_t cal_sig2(double sp[4], double E1, double E2)
{
	double tmp3 = sp[0]*E2+sp[1];
	double tmp4 = sp[2]*E1+sp[3];
	return sqrt(tmp3*tmp3*tmp4*tmp4/(tmp3*tmp3+tmp4*tmp4));
}

void MyMainFrame::DoFit()
{
	MyMainFrame::DoDraw();
	MyMainFrame::DoDraw2();
	
	TCanvas *fCanvas = fEcanvas2->GetCanvas();
	fCanvas->cd();
	
	// Lay cac tham so sigma (FWHM)
	stringstream ss_sig(bt_sig_par->GetText());
	double a_sig_par[4];
	ss_sig >> a_sig_par[0] >> a_sig_par[1] >> a_sig_par[2] >> a_sig_par[3];
	
	//Tao ham lam khop
	TF1 *ffit = new TF1("ffit",fit_function,bt_left2->GetNumber(), bt_right2->GetNumber(),3+max_peak_num*3);
	
	//Get Number of peak, Esum, bkg_deg
	int peak_num = 0;
	{
		double tmp;
		stringstream ss_pos(bt_peak_pos->GetText());
		while (ss_pos >> tmp) {peak_num++;}
	}	
	double E_sum = bt_Esum->GetNumber();
	int bkg_deg = int(bt_bkg_deg->GetNumber());
	
	//Set Initial parameter for fitting;
	for (int i=0; i<=bkg_deg; i++) ffit->SetParameter(i,0.01);
	for (int i=bkg_deg+1; i<3; i++) ffit->FixParameter(i,0.0);
	for (int i=3+peak_num*3; i<=2+max_peak_num*3; i++) ffit->FixParameter(i,1.0);
	for (int i=peak_num; i<=max_peak_num; i++) ffit->FixParameter(3+3*i,0.0);
	
	stringstream ss_pos1(bt_peak_pos->GetText());
	for (int i=1; i<=peak_num; i++)
	{
		double pos; ss_pos1 >> pos; ffit->SetParameter(4+3*(i-1),pos*1.0);
		ffit->FixParameter(5+3*(i-1),2.35*cal_sig(a_sig_par,pos*1.0,E_sum*1.0-pos*1.0));
	}
	
	hfit->Fit(ffit,"R");
	fCanvas->Update();
	
	cout<<"Testing ...... "<<peak_num<<endl;
}


MyMainFrame::~MyMainFrame()
{
   // Clean up used widgets: frames, buttons, layouthints
   Cleanup();
}

//*********************Program***************************

TH1F *loadspe(string input)
{
   double count; int channel;
   ifstream if1;
   TH1F *h = new TH1F("h", "Full Spectrum", 16384, 0., 16384.);
   for (int i = 1; i<=16384; i++) h->SetBinContent(i, 0.);
   if1.open(input.c_str());
   //int channel = 1;
   while (if1 >> channel >> count) {
      h->SetBinContent(channel, count);
      //channel++;
   }
   if1.close();
   return h;
}

void peakprogram()
{
   new MyMainFrame(gClient->GetRoot(), 200, 200);
}