// #include "LOpticsOpt.h"

class LOpticsOpt;

LOpticsOpt * opt;
UInt_t NPara = 0;
Double_t OldMatrixArray[10000]={-99};	//NPara
Bool_t free[10000]={kFALSE};//NPara

Test()
{
	gROOT->LoadMacro("LOpticsOpt.C+");

	opt = new LOpticsOpt();

	opt->LoadDataBase("db_L.vdc.dat.0");
	opt->Print();

	/////////////////////////////////////////////
	//array - matrix test
	opt->fCurrentMatrixElems = &(opt->fTMatrixElems);

	opt->Matrix2Array(OldMatrixArray,free);
	opt->Array2Matrix(OldMatrixArray);

	/////////////////////////////////////////////

// 	opt->Print();
	opt->SaveDataBase("db_L.vdc.dat.cp");

	///////////////////////////////////////////

	opt->LoadRawData("SieveData.f51"/*,10*/);

	opt->PrepareSieve();
	
	opt->VerifyMatrix_Sieve();
	
	cout <<opt->SumSquareDTh()<<endl;
	TCanvas * c1 =opt->CheckSieve();
// 	c1->Print(DestDataBase+".png","png");
	
}


TString SourceDataBase = "OptTh/db_L.vdc.dat.15.Seed";
TString DestDataBase 	= "OptTh/db_L.vdc.dat.16.test";
Bool_t OptTh = kTRUE;

// TString SourceDataBase = "OptPhi/db_L.vdc.dat.13.Seed";
// TString DestDataBase 	= "OptPhi/db_L.vdc.dat.14.MoreTerm.NoT310";
// Bool_t OptTh = kFALSE;

// TString SourceDataBase 	= "OptTh/db_L.vdc.dat.10.Test";
// TString DestDataBase 	= "OptTh/db_L.vdc.dat.10";
// Bool_t OptTh = kTRUE;


//____________________________________________________________________
void myfcn(Int_t &/*i*/, Double_t *, Double_t &f, Double_t *par, Int_t /*j*/)
{
   //minimisation function computing the sum of squares of residuals

	assert(opt);
	assert(opt->fCurrentMatrixElems);

	opt->Array2Matrix(par);

	if (OptTh)
		f= opt->SumSquareDTh();
	else
		f= opt->SumSquareDPhi();

	return;
}


void DoMinTP(void)
{
	///////////////////////////////////////////////////////////////////
// 	gROOT->LoadMacro("LOpticsOpt.C+");
	opt = new LOpticsOpt();
	
	///////////////////////////////////////////////////////////////////
	
	if (OptTh){
		cout << "Optimizing for Theta\n";
		opt->fCurrentMatrixElems = &(opt->fTMatrixElems);
		
	}
	else{
		cout << "Optimizing for Phi\n";
		opt->fCurrentMatrixElems = &(opt->fPMatrixElems);
	}
	
	///////////////////////////////////////////////////////////////////
	opt->LoadDataBase(SourceDataBase);
	NPara = opt->Matrix2Array(OldMatrixArray,free);
	opt->LoadRawData("SieveData.f51"/*,1000*/);
	opt->PrepareSieve();
	opt->VerifyMatrix_Sieve();


	///////////////////////////////////////////////////////////////////
	opt->Print();
	///////////////////////////////////////////////////////////////////
	TVirtualFitter::SetDefaultFitter("Minuit2");  //default is Minuit
	TVirtualFitter *fitter = TVirtualFitter::Fitter(NULL, NPara);
	fitter->SetFCN(myfcn);

	for(UInt_t i=0; i<NPara; i++)
	{
		Double_t absold = TMath::Abs(OldMatrixArray[i]);
		Double_t abslimit = absold>0?absold*10000:10000;
		
		fitter->SetParameter( i,Form("TMatrix%03d",i),OldMatrixArray[i],
							  absold>0?absold/10:0.1, -abslimit,abslimit);
// 		fitter->SetParameter(1,"asdf",0,0,0,0);

		if (!free[i]) fitter->FixParameter(i);
	}

	///////////////////////////////////////////////////////////////////
	fitter->Print();
	cout<<fitter->GetNumberFreeParameters()
			<<" Free  / "
			<<fitter->GetNumberTotalParameters()
			<<" Parameters\n";
	
	assert(opt->fNRawData>0);
	assert(NPara>0);
	assert(fitter->GetNumberFreeParameters()>0);
	assert(fitter->GetNumberTotalParameters()==NPara);
	
	
	///////////////////////////////////////////////////////////////////
	Double_t arglist[1] = {0};
	fitter->ExecuteCommand("MIGRAD", arglist, 0);

	///////////////////////////////////////////////////////////////////
	opt->Print();
	opt->SaveDataBase(DestDataBase);
	
	opt->SumSquareDTh();
	TCanvas * c1 =opt->CheckSieve();
	c1->Print(DestDataBase+".png","png");
// 	delete fitter;

// 	for (UInt_t i=0;i<100;i++) cout<<"\n";
	gSystem->Exec(Form("cp -vf %s %s.source"
			,SourceDataBase.Data(),DestDataBase.Data()));
	gSystem->Exec(Form("cp -vf LOpticsOptScript.log %s.log"
			,DestDataBase.Data()));
	
}

PlotSourceDataBase()
{
	PlotDataBase(SourceDataBase);
}

PlotDataBase(TString DatabaseFileName)
{
	

	opt = new LOpticsOpt();

	opt->LoadDataBase(DatabaseFileName);
	opt->Print();

	///////////////////////////////////////////

	opt->LoadRawData("SieveData.f51"/*,10*/);

	opt->PrepareSieve();
	opt->VerifyMatrix_Sieve();
	cout <<"SumSquareDTh = "<<opt->SumSquareDTh()<<endl;
	
	TCanvas * c1 =opt->CheckSieve();
	SaveCanvas(c1,DatabaseFileName+".sieve",0);
	
}

LOpticsOptScript()
{
	gROOT->LoadMacro("LOpticsOpt.C+");
	gStyle->SetOptStat(11);
  //turn off xy grids
	OKStyle->SetPadGridX(0);
	OKStyle->SetPadGridY(0);
	
	DoMinTP();
// 	PlotSourceDataBase();
// 	PlotDataBase("Production5/db_L.vdc.dat.Prod5.1.2GeV");
// 	PlotDataBase("OptPhi/db_L.vdc.dat.Huan.New");
}