#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;
   TGTextEntry *bt_spe_name;
   TGCheckButton *bt_logscale;
   TGNumberEntry *bt_left;
   TGNumberEntry *bt_right;
   TGNumberEntry *bt_left2;
   TGNumberEntry *bt_right2;
   TGCheckButton *bt_tscfit;
   TGCheckButton *bt_logscale2;
   TGCheckButton *bt_offset;
   TGCheckButton *bt_offset2;
   TGTextEntry *bt_peak_pos;
   TGTextEntry *bt_sig_par; // TSC n2g
   TGNumberEntry *bt_Esum; 
   TGNumberEntry *bt_bkg_deg; 
   TGNumberEntry *bt_ini_sig; //normal 
   TGTextEntry *bt_roi_name;
   TH1F *hspe;
   TH1F *hspe1;
   TH1F *hfit;
   TGLabel  *coordinate;
public:
   MyMainFrame(const TGWindow *p, UInt_t w, UInt_t h);
   virtual ~MyMainFrame();
   void DoDraw();
   void DoDraw2();
   void EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected);
   void AutoCal();
};

void MyMainFrame::EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected)
{
//  Writes the event status in the status bar parts

   const char *text0, *text1, *text3;
   char text2[50];
   //text0 = selected->GetTitle();
   //SetStatusText(text0,0);
   //text1 = selected->GetName();
   //SetStatusText(text1,1);
   //if (event == kKeyPress)
      // sprintf(text2, "%c", (char) px);
   // else
      // sprintf(text2, "%d,%d", px, py);
   //SetStatusText(text2,2);
   text3 = selected->GetObjectInfo(px,py);
   coordinate->SetText(text3);
   //SetStatusText(text3,3);
}

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_tscfit = new TGCheckButton(v11frame, "TSC_fit",1);
   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_tscfit, 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));
   
   bt_tscfit->SetState(1);
   
   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_sig2 = new TGLabel(v21frame, "Initial sigma");
	bt_ini_sig = new TGNumberEntry(v21frame,2.0);
	TGLabel *label_bkg = new TGLabel(v21frame, "Bkg poly degree (<=2)"); 
	bt_bkg_deg = new TGNumberEntry(v21frame,1,0,2);
	coordinate = new TGLabel(v21frame,"0,0");	
	
	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_sig2, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 0));
	v21frame->AddFrame(bt_ini_sig, 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(coordinate, 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);
	Int_t wid = fEcanvas2->GetCanvasWindowId();
    TCanvas *myc = new TCanvas("MyCanvas", 10,10,wid);
    fEcanvas2->AdoptCanvas(myc);
    myc->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)","MyMainFrame",this,
               "EventInfo(Int_t,Int_t,Int_t,TObject*)");
	
	v22frame->AddFrame(fEcanvas2, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY, 0, 0, 0, 0));
	
	//V23
	TGTextButton *bt_autocal = new TGTextButton(v23frame, "Auto_cal");
	bt_roi_name = new TGTextEntry(v23frame,"input.roi");
	v23frame->AddFrame(bt_roi_name, new TGLayoutHints(kLHintsExpandX, 4, 4, 4, 0));
	v23frame->AddFrame(bt_autocal, new TGLayoutHints(kLHintsExpandX, 4, 4, 10, 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()");
	
	bt_autocal->Connect("Clicked()", "MyMainFrame", this, "AutoCal()");
	
   // 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;
	if (bt_tscfit->IsOn())
	{
		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));
		}
	}
	else 
	{
		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->SetParameter(5+3*(i-1),2.35*bt_ini_sig->GetNumber());
		}
	}
	
	hfit->Fit(ffit,"R");
	fCanvas->Update();
	
	//Result Print
	Double_t par[1000];
	Double_t *upar;
	upar = new double[1000];
	ffit->GetParameters(par);
	upar = ffit->GetParErrors();
	
	printf("Fitting Result \n");
	printf("Fitting Range:  %i *********** %i \n",bt_left2->GetNumber(),bt_right2->GetNumber());
	printf("Number of peaks: %i \n",peak_num);
	printf("%10s%10s%15s%15s%8s%8s\n","POS","UPOS","AREA","UAREA","FWHM","UA(\%)");
	for (int i = 1; i<=peak_num; i++)
	{
		printf("%10.2f",par[4+3*(i-1)]+1.0);
		printf("%10.2f",upar[4+3*(i-1)]);
		printf("%15.2f",par[3+3*(i-1)]);
		printf("%15.2f",upar[3+3*(i-1)]);
		printf("%8.4f",par[5+3*(i-1)]);
		printf("%8.0f \n",100*upar[3+3*(i-1)]/par[3+3*(i-1)]);
	}
	//cout<<"Testing ...... "<<peak_num<<endl;
}

void MyMainFrame::AutoCal()
{
	ifstream if1;
	string s1 = bt_roi_name->GetText(); 
	if1.open(s1.c_str());
	string na_spename; if1 >> na_spename; 
	bt_spe_name->SetText(na_spename.c_str());
	if1 >> s1; 
	bt_tscfit->SetState(0);
	if (s1=="tsc") bt_tscfit->SetState(1);
	string s2=""; 
	for (int i=0; i<4; i++) {if1 >> s1; if (i<3) s2 = s2+s1+" "; if (i==3) s2=s2+s1;} 
	bt_sig_par->SetText(s2.c_str());
	double d1; 
	if1 >> d1; 
	bt_Esum->SetNumber(d1);
	int i1, i2, i3; 
	if1 >> i1 >> i2 >> i3;
	bt_left2->SetNumber(i1);
	bt_right2->SetNumber(i2);
	bt_peak_pos->SetText(""); s2 ="";
	for (int i=0; i<i3; i++)
	{
		if1 >> s1;
		if (i<i3-1) s2 = s2 + s1 + " ";
		if (i==i3-1) s2 = s2 + s1; 
	}
	bt_peak_pos->SetText(s2.c_str());
	//MyMainFrame::DoFit();
	//cout<< na_spename<<endl;
	if1.close();
	
}


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);
}