// #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 	= "DBchange/db_L.vdc.dat.010";
TString DestDataBase 	= "DBchange/db_L.vdc.dat.012";
TString RawDataFile     = "Dp/Dp_LH2.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->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(RawDataFile);
	opt->PrepareDp();
	//opt->VerifyMatrix_Sieve();
	opt->VerifyMatrix_Dp();

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

	opt->fArbitaryDpKinShift[0] = 8.930709e-05;
	opt->fArbitaryDpKinShift[1] = 9.338335e-05;
	opt->fArbitaryDpKinShift[2] = -3.086960e-06;
	opt->fArbitaryDpKinShift[3] = 1.429621e-04;
	opt->fArbitaryDpKinShift[4] = 1.043272e-04;
	
	///////////////////////////////////////////////////////////////////
	opt->Print();
	///////////////////////////////////////////////////////////////////
	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->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 <<"SumSquareDp = "<<opt->SumSquareDp()<<endl;
	
	TCanvas * c1 =NULL;

	c1 =opt->CheckDpVSAngle();
	SaveCanvas(c1,DatabaseFileName+".VSAngle");
	
	c1 =opt->CheckDpVSCutID();
	SaveCanvas(c1,DatabaseFileName+".VSID");
	
	c1 =opt->CheckDp();
	SaveCanvas(c1,DatabaseFileName+".dp",0);
}

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("OptDp/db_L.vdc.dat.17.MidPlnSym.RevRadLoss.ReOptimized");
// 	PlotDataBase("Production5/db_L.vdc.dat.Prod5.1.2GeV");
}