#include"TROOT.h"
#include"TMath.h"
#include<iostream>
#include"TF1.h"
#define N 580
using namespace std;
Double_t z1[N];
Double_t u[N];
Double_t s[N];
Double_t mod(Double_t* z,Double_t* par)
{
return sqrt(1./((4.173558096*1.e-5/(par[0]*par[0]))*(1+z[0])*(1+z[0])*(1+z[0])*(1+z[0])+par[1]*(1+z[0])*(1+z[0])*(1+z[0])+(1-(4.173558096*1.e-5/(par[0]*par[0]))-par[1])*pow(1+z[0],3.*(1.+par[2]))));
}
Double_t magnetude(Double_t* z,Double_t* par)
{ 
TF1 f("f",mod,0,1.75,3);
f.SetParameters(par[0],par[1],par[2]);
//return 5.*TMath::Log10(((1.+z[0])*(f.Integral(0.,z[0])))/par[0])+42.38;
return 5.*TMath::Log10(((1.+z[0])*(3000./par[0])*(f.Integral(0.,z[0]))))+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;
   Double_t err=0;

   for (i=0;i<nbins; i++) 
   {
     delta  = double(u[i]-magnetude(&z1[i],par))/s[i];
     //err=1./(s[i]*s[i]+0.15*0.15+0.01*0.01);
     chisq += delta*delta;
   }
   f = chisq;
}
//_________________________________________________________________________________________
void BR()
{
double _z,_u,_s;
ifstream file ("union2.1");
for(int i=0;i<N;i++)
{
file>>_z>>_u>>_s;
z1[i]=_z; u[i]=_u; s[i]=_s;
}
////////////
   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 ,4,ierflg);
   arglist[0] = 2;
   gMinuit->mnexcm("SET STR", arglist ,4,ierflg);
   
   //arglist[0] = 0.1;


   //gMinuit->mnexcm("SET EPS", arglist, 1, ierflg); 
// Set starting values and step sizes for parameters
   gMinuit->mnparm(0, "h   ", 0.7,0.01, 0,0,ierflg);
   gMinuit->mnparm(1, "omega_m", 0.3,0.01, 0,0,ierflg);
   gMinuit->mnparm(2, "wd", -1.02,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.001;//tolerance
   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
   gMinuit->mnexcm("MINOS", arglist ,2,ierflg);
   //gMinuit->mnimpr();
// Print results
  Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);   
gMinuit->mnprin(4,amin);

}