// #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 	= "OptDp/db_L.vdc.dat.6.Corrph";
TString DestDataBase 	= "OptDp/db_L.vdc.dat.6.Corrph.OffSet";

//____________________________________________________________________
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->SumSquareDp();

	return;
}

void DoMin(void)
{
	///////////////////////////////////////////////////////////////////
// 	gROOT->LoadMacro("LOpticsOpt.C+");
	opt = new LOpticsOpt();
	
	///////////////////////////////////////////////////////////////////
	
	cout << "Optimizing for dp\n";
	opt->fCurrentMatrixElems = &(opt->fDMatrixElems);
	
	///////////////////////////////////////////////////////////////////
	opt->LoadDataBase(SourceDataBase);
	NPara = opt->Matrix2Array(OldMatrixArray,free);
	opt->LoadRawData("DpData.f51"/*,1000*/);
	opt->PrepareDp();
	opt->VerifyMatrix_Sieve();
	opt->VerifyMatrix_Dp();

	///////////////////////////////////////////////////////////////////
	//compensate bias due to dp event selections

	opt->fArbitaryDpKinShift[0] = 1.144078e-04;
	opt->fArbitaryDpKinShift[1] = 8.856055e-05;
	opt->fArbitaryDpKinShift[2] = -2.194187e-05;
	opt->fArbitaryDpKinShift[3] = 1.245723e-04;
	opt->fArbitaryDpKinShift[4] = 1.130130e-04;
	
	///////////////////////////////////////////////////////////////////
	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->SumSquareDp();
	
	TCanvas * c1 =opt->CheckDp();
	c1->Print(DestDataBase+".png","png");
	
	TCanvas * c2 =opt->CheckDpVSAngle();
	c2->Print(DestDataBase+".VSAngle.png","png");
	
	TCanvas * c3 =opt->CheckDpVSCutID();
	c3->Print(DestDataBase+".VSID.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()));
}

PlotDataBase(TString DatabaseFileName, Bool_t WithAngle=kFALSE)
{
	opt = new LOpticsOpt();

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

	///////////////////////////////////////////
	opt->LoadRawData("DpData.f51"/*,1*/);

	opt->PrepareDp();
	opt->VerifyMatrix_Sieve();
	opt->VerifyMatrix_Dp();
	cout <<"SumSquareDTgY = "<<opt->SumSquareDp()<<endl;
	
	TCanvas * c1 =NULL;

	c1 =opt->CheckDpVSAngle();
	c1->Print(DatabaseFileName+".VSAngle.png","png");
	
	c1 =opt->CheckDpVSCutID();
	c1->Print(DatabaseFileName+".VSID.png","png");
	
	c1 =opt->CheckDp();
	c1->Print(DatabaseFileName+".png","png");

	
	
}

PlotSourceDataBase( Bool_t WithAngle=kFALSE)
{
	PlotDataBase(SourceDataBase,WithAngle);
}

PlotDestDataBase( Bool_t WithAngle=kFALSE)
{
	PlotDataBase(DestDataBase,WithAngle);
}

Test()
{
	
	opt = new LOpticsOpt();
	opt->LoadDataBase("OptDp/db_L.vdc.dat.0");
	opt->Print();
	
	opt->LoadRawData("DpData.f51"/*,10*/);

	opt->PrepareDp();
	opt->VerifyMatrix_Sieve();
	opt->VerifyMatrix_Dp();
// 	cout <<"SumSquareDTgY = "<<opt->SumSquareDTgY()<<endl;
	
	TCanvas * c1 =opt->CheckDp();
	c1->Print(DatabaseFileName+".png","png");
}

LOpticsOptScriptDp()
{//choose the entracne function here
	gROOT->LoadMacro("LOpticsOpt.C+");
	gStyle->SetOptStat(11);
	
// 	DoMin();
// 	Test();
// 	PlotSourceDataBase();
// 	PlotDestDataBase(kTRUE);
	PlotDataBase("Production4/db_L.vdc.dat.Prod4");
// 	PlotDataBase("OptDp/db_L.vdc.dat.Huan.New");
}