#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 LRmod(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])+(9./4.)*par[2]*(TMath::Log(1+z[0]))*(TMath::Log(1+z[0]))-3.*sqrt(par[2]*(1-(4.173558096*1.e-5/(par[0]*par[0]))-par[1]))*TMath::Log(1+z[0])+(1-((4.173558096*1.e-5)/(par[0]*par[0]))-par[1])));
}

Double_t magnetude(Double_t* z,Double_t* par)
{ 
TF1 f("f",LRmod,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;
  
   for (i=0;i<nbins; i++) 
   {
     delta  = (u[i]-magnetude(&z1[i],par))/s[i];

     chisq += delta*delta;
   }
   f = chisq;
}
//_________________________________________________________________________________________
void LR()
{
Double_t _z=0,_u=0,_s=0;
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.01;

   //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, "om*", 1.e-4,1.e-5, 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("HESSE", 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);//print function value, estimated distance to minimum, parameter errors, step sizes 

}