//
//Polynomial fit
//
/*
For 2 par fit : p0 = 3.2851e-10
                p1 = 9.4419e-10, matched well.

For 6 par fit, result sould be
   p0 :  3.275E-10
   p1 :  2.735E-09
   p2 :  1.760E-10

   p3 :  9.494E-10
   p4 : -1.700E-10
   p5 :  5.803E-09
*/

//

const Int_t n = 1522;

Double_t x[n];
Double_t y[n];

Double_t xerr[n];
Double_t yerr[n];

//______________________________________________________________________________
void fit2pv1()
{

	Double_t xin, yin, xein, yein;

	//
	ifstream in;
	in.open("sample.txt");

	Int_t n = 0;
	while (1) {

		in >> xin >> yin >> xein >> yein;

		x[n] = xin*1.00e-03; // Unit : xin,yin   : [mm]
		y[n] = yin*1.00e-03; //        xein,yein : [um]

		xerr[n] = xein*1.00e-06;
		yerr[n] = yein*1.00e-06;

		//
		if (!in.good()) break;
		if (n%100==0) printf("%d-th : x=%8f, y=%8f, xerr=%8f, yerr=%8f\n",n,xin,yin,xein,yein);
		n++;
	}
	in.close();
	cout << endl;

	// TMinuit
	TMinuit *gMinuit = new TMinuit(2);  //initialize TMinuit
	gMinuit->SetFCN(fcn);

	Double_t arglist[2];
	Int_t ierflg = 0;

	arglist[0] = 1;
	gMinuit->mnexcm("SET ERR",arglist,1,ierflg);


	// Set starting values and step sizes for parameters
	Double_t vstart[2] = {1.00e-10, 1.00e-10};
	Double_t step[2]   = {1.00e-12, 1.00e-12};
	
	gMinuit->mnparm(0, "p0", vstart[0], step[0], 0,0,ierflg);
	gMinuit->mnparm(1, "p1", vstart[1], step[1], 0,0,ierflg);

	// Now ready for minimization step
	arglist[0] = 500;
	arglist[1] = 1.;

	cout << "-- MIGRAD --" << endl;
	gMinuit->mnexcm("MIGRAD",arglist,1,ierflg);

	// Print results
	Double_t amin,edm,errdef;
	Int_t nvpar,nparx,icstat;

	gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
	gMinuit->mnprin(3,amin);
	cout << endl;

	// Check MINOS
	//cout << "-- MINOS --" << endl;
	//gMinuit->mnexcm("MINOS",arglist,1,ierflg);

}


//______________________________________________________________________________
Double_t funcX(Double_t x, Double_t y, Double_t *par)
{

  Double_t value = par[0]*TMath::Power(x,0);
  return value;
}

//______________________________________________________________________________
Double_t funcY(Double_t x, Double_t y, Double_t *par)
{

  Double_t value = par[1]*TMath::Power(y,0);
  return value;
}

//______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{

	Int_t i;
	Double_t chisq = 0.;

	Double_t dx, dy;

	for (i=0; i < n; i++) {

		dx = funcX(x[i]+xerr[i],y[i]+yerr[i],par) - xerr[i];
		dy = funcY(x[i]+xerr[i],y[i]+yerr[i],par) - yerr[i];

		chisq += dx*dx + dy*dy;
	}
	f = chisq;

}