// #include "LOpticsOpt.h"

class LOpticsOpt;

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

Int_t Flag=1;//1=optimize DTgY,2=optimize DPhi, 3=optimize DTh


TString SourceDataBase 	= "DBchange/db_L.vdc.dat.006";
TString DestDataBase 	= "DBchange/db_L.vdc.dat.010";
if(Flag==2 || Flag==3) TString RawDataFile     = "Sieve/Sieve_1893.f51";
if(Flag==1) TString RawDataFile     = "Vertex/Vertex_2013.f51";


//____________________________________________________________________
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(Flag==2) f=opt->SumSquareDPhi();
	if(Flag==3) f=opt->SumSquareDTh();

	if(Flag==1) f= opt->SumSquareDTgY();
 	if(Flag==4) f= opt->SumSquareDTgYAverFoils();

	return;
}

void DoMin(void)
{
  ///////////////////////////////////////////////////////////////////
  // 	gROOT->LoadMacro("LOpticsOpt.C+");
  opt = new LOpticsOpt();
  
  ///////////////////////////////////////////////////////////////////
  
  if(Flag==1){
    cout << "Optimizing for tg_y\n";
    opt->fCurrentMatrixElems = &(opt->fYMatrixElems);
  }
  if(Flag==2){
    cout << "Optimizing for tg_phi\n";
    opt->fCurrentMatrixElems = &(opt->fPMatrixElems);
  }
  if(Flag==3){
    cout << "Optimizing for tg_th\n";
    opt->fCurrentMatrixElems = &(opt->fTMatrixElems);
  }
	///////////////////////////////////////////////////////////////////



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

	opt->LoadDataBase(SourceDataBase);
	opt->Print();
	NPara = opt->Matrix2Array(OldMatrixArray,free);
	opt->LoadRawData(RawDataFile/*,1000*/);
// 	opt->LoadRawData("VertexDataRasterBeamWithPointing.f51"/*,1000*/);
	//opt->PrepareVertex();
	opt->VerifyMatrix_Vertex();

	opt->PrepareSieve();
	opt->VerifyMatrix_Sieve();
	//opt->VerifyMatrix_Vertex();

	///////////////////////////////////////////////////////////////////
	TVirtualFitter::SetDefaultFitter("Minuit");  //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->SumSquareDTgY();
	opt->SumSquareDTh();
	//TCanvas * c1 =opt->CheckVertex();
	TCanvas * c1 =opt->CheckSieve();
	c1->Print(DestDataBase+".png","png");
// 	delete fitter;

	opt->SaveDataBuffer(DestDataBase+".root");

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

PlotDestDataBase()
{
	PlotDataBase(DestDataBase);
}

PlotDataBase(TString DatabaseFileName)
{
	opt = new LOpticsOpt();

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

	///////////////////////////////////////////
	opt->LoadRawData(RawDataFile/*,1*/);
	//opt->PrepareVertex();
	opt->PrepareSieve();
	opt->VerifyMatrix_Sieve();
	//opt->VerifyMatrix_Vertex();
	//cout <<"SumSquareDTgY = "<<opt->SumSquareDTgY()<<endl;
	cout <<"SumSquareDTh = "<<opt->SumSquareDTh()<<endl;
	
	//TCanvas * c1 =opt->CheckVertex();
// 	c1->Print(DatabaseFileName+".Vertex.png","png");
// 	c1->Print(DatabaseFileName+".Vertex.C","cxx");
	TCanvas * c1 =opt->CheckSieve();
	c1->Print(DatabaseFileName+".Sieve.png","png");
	c1->Print(DatabaseFileName+".Sieve.C","cxx");
	
}

LOpticsOptScriptVertex()
{
	gROOT->LoadMacro("LOpticsOpt.C+");
	gStyle->SetOptStat(11);
// 	
 	DoMin();
	//PlotSourceDataBase();
	//PlotDestDataBase();
// 	PlotDataBase("OptVZ/db_L.vdc.dat.6");
//	PlotDataBase("Production5/db_L.vdc.dat.Prod5.2.4GeV");
}