// #include "LOpticsOpt.h"

class LOpticsOpt;

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

// TString SourceDataBase 	= "OptVZ/db_L.vdc.dat.11.Seed";
// TString DestDataBase 	= "OptVZ/db_L.vdc.dat.12.OffSet3";
// TString RawDataFile = "VertexDataRasterBeamWithPointing.f51";


TString SourceDataBase 	= "OptVZ2505/db_L.vdc.dat.15.AddY400";
TString DestDataBase 	= "OptVZ2505/db_L.vdc.dat.15.AddY400.OffSet2";
TString RawDataFile = "VertexData2505.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);
	
	f= opt->SumSquareDTgY();
// 	f= opt->SumSquareDTgYAverFoils();

	return;
}

void DoMin(void)
{
	///////////////////////////////////////////////////////////////////
// 	gROOT->LoadMacro("LOpticsOpt.C+");
	opt = new LOpticsOpt();
	
	///////////////////////////////////////////////////////////////////
	
	cout << "Optimizing for tg_y\n";
	opt->fCurrentMatrixElems = &(opt->fYMatrixElems);
	
	///////////////////////////////////////////////////////////////////
	opt->fArbitaryVertexShift[0] = 1.331942e-03;
	opt->fArbitaryVertexShift[1] = -6.610745e-04;
	opt->fArbitaryVertexShift[2] = -6.046560e-04;
	opt->fArbitaryVertexShift[3] = -2.443946e-04;
	opt->fArbitaryVertexShift[4] = -2.505168e-04;
	opt->fArbitaryVertexShift[5] = -6.914623e-04;
	opt->fArbitaryVertexShift[6] = -8.919399e-04;

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

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

	///////////////////////////////////////////////////////////////////
	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->SumSquareDTgY();
	TCanvas * c1 =opt->CheckVertex();
	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);
}
PlotDestDataBase()
{
	PlotDataBase(DestDataBase);
}

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

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

	///////////////////////////////////////////
	opt->LoadRawData(RawDataFile/*,1*/);

	opt->PrepareVertex();
	opt->VerifyMatrix_Sieve();
	opt->VerifyMatrix_Vertex();
	cout <<"SumSquareDTgY = "<<opt->SumSquareDTgY()<<endl;
	
	TCanvas * c1 =opt->CheckVertex();
	c1->Print(DatabaseFileName+".png","png");
	
}

LOpticsOptScriptVertex()
{
	gROOT->LoadMacro("LOpticsOpt.C+");
	gStyle->SetOptStat(11);
	
	DoMin();
// 	PlotSourceDataBase(); 
// 	PlotDataBase("OptVZ/db_L.vdc.dat.6");
// 	PlotDataBase("Production3/db_L.vdc.dat.Prod3");