#include"TROOT.h"
#include <iomanip>
#include <iostream>
#define N 580
using namespace std;
Double_t z1[N];
Double_t u[N];
Double_t s[N];
Double_t Fsab(Double_t* z,Double_t* par)
{

return 1./sqrt(par[0]*(1.+z[0])*(1.+z[0])*(1.+z[0])+(1.-par[0])*(1.+(par[1]/(1.-par[0]))*log(1.+z[0])));

}

Double_t magnetude(Double_t* z,Double_t* par)
{ 
TF1 f("f",Fsab,0,1.75,2);
f.SetParameters(par[0],par[1]);
//return 5.*TMath::Log10(((1.+z[0])*(f.Integral(0.,z[0])))/par[2])+42.38;
//return 43.3+5*TMath::Log10((1+z[0])*(f.Integral(0,z[0])));

return 5.*TMath::Log10(((1.+z[0])*3000.*(f.Integral(0.,z[0])))/par[2])+25.;
}


//__________________________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   const Int_t nbins = N;
   Int_t i;

//calculate chisquare
   Double_t chisq = 0.;
   Double_t delta=0.;
   for (i=0;i<nbins; i++) 
   {
    //if(s[i]!=0)
     delta  = (u[i]-magnetude(&z1[i],par))/s[i];
     chisq += delta*delta;
   }
}
//_________________________________________________________________________________________
void sab()
{
double _z,_u,_s;
ifstream file ("union2.1");
for(int i=0;i<N;i++)
{
file>>_z>>_u>>_s;
//cout.precision(15);
z1[i]=_z; u[i]=_u; s[i]=_s;
//cout<<z1[i]<<"::"<<u[i]<<"::"<<s[i]<<endl;
}

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

  TMinuit *gMinuit = new TMinuit(5);  //initialize TMinuit with a maximum of 5 params
   gMinuit->SetFCN(fcn);

   Double_t arglist[10];
   Int_t ierflg = 0;

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

   //gMinuit->mnexcm("SET EPS", arglist, 1, ierflg);

// Set starting values and step sizes for parameters
   //gMinuit->FixParameter(1);
   
   gMinuit->mnparm(0, "omega_m", 0.300, 0.01, 0,0,ierflg);
   gMinuit->mnparm(1, "A    ",1.e-6, 1.e-6,0,0,ierflg);
   gMinuit->mnparm(2, "h      ",0.7, 0.01,0,0,ierflg);
   
// Now ready for minimization step
//Call MIGRAD with 500 iterations maximum
   arglist[0] = 5000; //max number of iterations
   arglist[1] = 0.01; //tolerance
   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
   gMinuit->mnexcm("MINOS", arglist ,2,ierflg);
// Print results
  Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
gMinuit->mnprin(4,amin);

/*gMinuit->Command("SCAn 1 "); 
 
TGraph *gr = (TGraph*)gMinuit->GetPlot(); 
gr->SetMarkerStyle(21); 
gr->Draw("alp");*/

/*gMinuit->SetErrorDef(2.3);
TGraph *gr2 = (TGraph*)gMinuit->Contour(40,1,0);
gr2->SetFillColor(42);
gr2->GetXaxis()->SetTitle("A");
gr2->GetXaxis()->SetLimits(0.0,0.6); 
gr2->GetYaxis()->SetTitle("#Omega_{m}");
gr2->GetYaxis()->SetLimits(0,0.7);
gr2->Draw("alf");*/
}