/*
	Analyse stars around burst and output their positions/scale. From these stars with stable position can be chosen
	for burst analysis in anal_burst
*/

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "TH2D.h"
#include "TProfile2D.h"
#include "TGraph2D.h"
#include "TTree.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TSQLServer.h"
#include "TSQLResult.h"
#include "TSQLRow.h"
#include "TF1.h"
#include "TF2.h"
#include "TStyle.h"
#include "TTimer.h"
#include "TApplication.h"
#include "TVectorD.h"
#include "TEventList.h"
#include "TNtupleD.h"
#include "TMath.h"
#include "Math/MinimizerOptions.h"
#include "zernike_fit_area_speedup.h"
#include <vector>

using namespace std;

			TString answ;
		
			TTimer *timer =	new TTimer("gSystem->ProcessEvents();", 50, kFALSE);		

int main(int argc, char **argv)
{
	cerr.precision(12);
	cout.precision(11);

	Int_t old_argc=argc;
	
	TApplication theApp("fit_stars", &argc, argv);	
	
	gStyle->SetPalette(1,0);
	gStyle->SetOptStat(1111);
	gStyle->SetNumberContours(99);
	Float_t posx, posy, ccdx, ccdy, scale, mag, ra, dec, chi2, r, offset, axis;

	bool draw=true;

	int range=1;

	TCanvas *br_canv = new TCanvas("br_canv", "br_canv", 540, 500);
	TCanvas *c1 = new TCanvas("c1", "c1", 1050, 540);	
	TCanvas *c = new TCanvas("c", "c", 1050, 540);
	c->Divide(2,1);
	c1->Divide(2,1);

	
	// Load PRF
	TFile *f1 = new TFile("par_histos.root");
	prf_histo = (TH2D*)f1->Get("PRF");
	prf_histo->SetDirectory(0);
	f1->Close();

	TDatime *tm = new TDatime();

	// Initialize fitting	
	TF2 *mpsf = init_mpsf(-0.3926991, 1400, 6.5);
	mpsf->SetNpx(41);
	mpsf->SetNpy(41);
	mpsf->SetRange(-20.5, -20.5, 20.5, 20.5);
	
	ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLMultiFit", "");
//		ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Genetic2c", "");
//		ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "");
//		ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Fumili", "");
	ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(0);
	ROOT::Math::MinimizerOptions::SetDefaultStrategy(1);
//	ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.001);
//	ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.000001);
//	ROOT::Math::MinimizerOptions::SetDefaultPrecision(0.000000001);
//	ROOT::Math::MinimizerOptions::SetDefaultPrecision(1);
//		ROOT::Math::MinimizerOptions::SetDefaultMaxIterations(25);
	ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.001);
//	ROOT::Math::MinimizerOptions::SetDefaultTolerance(10);

	TH1D *br_histo=NULL;

	ccdx = 1363;
	ccdy = 925;

	axis = TMath::ATan2((ccdy-1031),(ccdx-1126))-TMath::Pi()/2;
	if(axis<0) axis+=TMath::TwoPi();
	r = sqrt((ccdx-1126)*(ccdx-1126)+(ccdy-1031)*(ccdy-1031));
	Int_t dax = axis*100000;
	axis = dax/100000.;

	// Store positions for this star

	range=2;
	scale_fit_range=range;
	
	mpsf->FixParameter(0, r);
	mpsf->FixParameter(3, axis);
	mpsf->FixParameter(6, 2.);
	mpsf->FixParameter(7, 1.);
	mpsf->FixParameter(8, 1);
	mpsf->FixParameter(4, 1);

	br_histo = new TH1D("br_histo", "br_histo", 6721, -120.5, 6600.5);

	// Fill the star histogram
	TFile *fhisto = new TFile("imghisto.root", "READ");
	img_histo = (TH2D*) fhisto->Get("img_histo");
	nbx = img_histo->GetNbinsX();
	nby = img_histo->GetNbinsY();			
	nbxo2=nbx/2.;
	nbyo2=nby/2.;
	img_histo->SetStats(false);

	for(int i=1; i<=img_histo->GetNbinsX(); i++) for(int j=1; j<=img_histo->GetNbinsX(); j++) br_histo->Fill(img_histo->GetBinContent(i,j));

	Long64_t med_cnt=0;
	Int_t med[41*41];

	// Calculate median in the +-3sigma range
	for(int i=1; i<=br_histo->GetEntries(); i++)
	{
		for(int j=1; j<=br_histo->GetBinContent(i); j++)
		{
			Int_t br_val = br_histo->GetBinCenter(i);
			if(br_val>(br_histo->GetMean()+br_histo->GetRMS()*2) || br_val<(br_histo->GetMean()-br_histo->GetRMS()*2)) continue;
			med[med_cnt] = br_val;
			med_cnt++;
		}
	}
	
	Int_t br_median = TMath::Median(med_cnt, med);
	
	if(draw) cout << "median: " << br_median << endl;			

	if(draw)
	{
		// Draw the star
		TString tit = TString::Format("%.0f %.0f", ccdx, ccdy);
		c->cd(1);
		gPad->SetRightMargin(0.14);			
		img_histo->DrawCopy("colz");
		gPad->Modified(); gPad->Update();
	}

	Double_t h_max = img_histo->GetMaximum();
	Double_t h_min = img_histo->GetMinimum();
	
	// Get image brightness histogram
	if(draw)
	{
		br_canv->cd();
		br_histo->Draw();
		gPad->Modified(); gPad->Update();
	}

	mpsf->FixParameter(4, 1);
	mpsf->FixParameter(5, br_median);


	// *** Fit the main star ***

	TFitResultPtr fptr;

	if(draw)
	{
		c->cd(2);
		mpsf->Draw("colz");
		gPad->Modified();
		gPad->Update();			
		timer->TurnOn();
		timer->Reset();

		cout << "Press a key and enter" << endl;
		cin >> answ;

		timer->TurnOff();				
	}

	Double_t init_x=0, init_y=0;

	// Not sure why assymetric limits work better for lm
	init_x = -0.15869591879;
	init_y = -0.3446373296;

	mpsf->SetParameter(1,init_x);
	mpsf->SetParameter(2,init_y);
	mpsf->SetParLimits(1, init_x-1.5, init_x+1.5);
	mpsf->SetParLimits(2, init_y-1.5, init_y+1.5);
//				mpsf->SetParError(1, 100);
//				mpsf->SetParError(2, 100);
				
	if(draw)
	{
		c->cd(2);
		mpsf->Draw("colz");
		gPad->Modified();
		gPad->Update();			
		timer->TurnOn();
		timer->Reset();

		cout << "Press a key and enter" << endl;
		cin >> answ;

		timer->TurnOff();				
		
	}

	cout << "Initial params " << endl;
	mpsf->Print();
	
	fptr = img_histo->Fit(mpsf, "RNSV");
}
