#include"TROOT.h"
#include"TMath.h"
#include<iostream>
#include"TF1.h"
#include"TMatrixD.h"
#define N 580
using namespace std;
Double_t z1[N];
Double_t u[N];
Double_t s[N];
//TH1F* h= new TH1F("hx","Distrubution X2",100,500,590);
Double_t mod(Double_t* z,Double_t* par)
{
return sqrt(1./(par[0]*(1+z[0])*(1+z[0])*(1+z[0])+(1-par[0])*pow(1+z[0],3.*(1.+par[1]))));
//return sqrt(1./(par[0]*(1+z[0])*(1+z[0])*(1+z[0])*(1-pow(1+z[0],3.*par[1]))+pow(1+z[0],3.*(1.+par[1]))));

}
Double_t magnetude(Double_t* z,Double_t* par)
{ 
TF1 f("f",mod,0.015,1.414,2);
f.SetParameters(par[0],par[1]);
//return 5.*TMath::Log10(((1.+z[0])*(f.Integral(0.,z[0])))/par[0])+42.38;
//cout<<"integral="<<f.Integral(0.,z[0]))<<endl;
return 5.*TMath::Log10(((3000.*(1.+z[0]))/par[2])*(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  = (u[i]-magnetude(&z1[i],par))/s[i];
     chisq += delta*delta;
   }
   f = chisq;
}
//_________________________________________________________________________________________
void variable()
{
double var[580];
ifstream file ("mat21.txt");
ofstream res("fcn.txt");
//TMatrixD mat(580,580);
for(int i=0;i<N;i++)
{
file>>var[0]>>var[1]>>var[2]>>var[3]>>var[4]>>var[5]>>var[6]>>var[7]>>var[8]>>var[9]>>var[10]>>var[11]>>var[12]>>var[13]>>var[14]>>var[15]>>var[16]>>var[17]>>var[18]>>var[19]>>var[20]>>var[21]>>var[22]>>var[23]>>var[24]>>var[25]>>var[26]>>var[27]>>var[28]>>      var[29]>>var[30]>>var[31]>>var[32]>>var[33]>>var[34]>>var[35]>>var[36]>>var[37]>>var[38]>>var[39]>>var[40]>>var[41]>>var[42]>>var[43]>>var[44]>>var[45]>>var[46]>>var[47]>>var[48]>>var[49]>>var[50]>>var[51]>>var[52]>>var[53]>>var[54]>>var[55]>>var[56]>>var[57]>>var[58]>>var[59]>>var[60]>>var[61]>>var[62]>>var[63]>>var[64]>>var[65]>>var[66]>>var[67]>>var[68]>>var[69]>>var[70]>>var[71]>>var[72]>>var[73]>>var[74]>>var[75]>>var[76]>>var[77]>>var[78]>>var[79]>>var[80]>>var[81]>>var[82]>>var[83]>>var[84]>>var[85]>>var[86]>>var[87]>>var[88]>>var[89]>>var[90]>>var[91]>>var[92]>>var[93]>>var[94]>>var[95]>>var[96]>>var[97]>>var[98]>>var[99]>>var[100]>>
var[101]>>var[102]>>var[103]>>var[104]>>var[105]>>var[106]>>var[107]>>var[108]>>var[109]>>var[110]>>var[111]>>var[112]>>var[113]>>var[114]>>var[115]>>var[116]>>var[117]>>var[118]>>var[119]>>var[120]>>var[121]>>var[122]>>var[123]>>var[124]>>var[125]>>var[126]>>var[127]>>var[128]>>var[129]>>var[130]>>var[131]>>var[132]>>var[133]>>var[134]>>var[135]>>var[136]>>var[137]>>var[138]>>var[139]>>var[140]>>var[141]>>var[142]>>var[143]>>var[144]>>var[145]>>var[146]>>var[147]>>var[148]>>var[149]>>var[150]>>var[151]>>var[152]>>var[153]>>var[154]>>var[155]>>var[156]>>var[157]>>var[158]>>var[159]>>var[160]>>var[161]>>var[162]>>var[163]>>var[164]>>var[165]>>var[166]>>var[167]>>var[168]>>var[169]>>var[170]>>var[171]>>var[172]>>var[173]>>var[174]>>var[175]>>var[176]>>var[177]>>var[178]>>var[179]>>var[180]>>var[181]>>var[182]>>var[183]>>var[184]>>var[185]>>var[186]>>var[187]>>var[188]>>var[189]>>var[190]>>var[191]>>var[192]>>var[193]>>var[194]>>var[195]>>var[196]>>var[197]>>var[198]>>var[199]>>var[200]>>
var[201]>>var[202]>>var[203]>>var[204]>>var[205]>>var[206]>>var[207]>>var[208]>>var[209]>>var[210]>>var[211]>>var[212]>>var[213]>>var[214]>>var[215]>>var[216]>>var[217]>>var[218]>>var[219]>>var[220]>>var[221]>>var[222]>>var[223]>>var[224]>>var[225]>>var[226]>>var[227]>>var[228]>>var[229]>>var[230]>>var[231]>>var[232]>>var[233]>>var[234]>>var[235]>>var[236]>>var[237]>>var[238]>>var[239]>>var[240]>>var[241]>>var[242]>>var[243]>>var[244]>>var[245]>>var[246]>>var[247]>>var[248]>>var[249]>>var[250]>>var[251]>>var[252]>>var[253]>>var[254]>>var[255]>>var[256]>>var[257]>>var[258]>>var[259]>>var[260]>>var[261]>>var[262]>>var[263]>>var[264]>>var[265]>>var[266]>>var[267]>>var[268]>>var[269]>>var[270]>>var[271]>>var[272]>>var[273]>>var[274]>>var[275]>>var[276]>>var[277]>>var[278]>>var[279]>>var[280]>>var[281]>>var[282]>>var[283]>>var[284]>>var[285]>>var[286]>>var[287]>>var[288]>>var[289]>>var[290]>>var[291]>>var[292]>>var[293]>>var[294]>>var[295]>>var[296]>>var[297]>>var[298]>>var[299]>>var[300]>>
var[301]>>var[302]>>var[303]>>var[304]>>var[305]>>var[306]>>var[307]>>var[308]>>var[309]>>var[310]>>var[311]>>var[312]>>var[313]>>var[314]>>var[315]>>var[316]>>var[317]>>var[318]>>var[319]>>var[320]>>var[321]>>var[322]>>var[323]>>var[324]>>var[325]>>var[326]>>var[327]>>var[328]>>var[329]>>var[330]>>var[331]>>var[332]>>var[333]>>var[334]>>var[335]>>var[336]>>var[337]>>var[338]>>var[339]>>var[340]>>var[341]>>var[342]>>var[343]>>var[344]>>var[345]>>var[346]>>var[347]>>var[348]>>var[349]>>var[350]>>var[351]>>var[352]>>var[353]>>var[354]>>var[355]>>var[356]>>var[357]>>var[358]>>var[359]>>var[360]>>var[361]>>var[362]>>var[363]>>var[364]>>var[365]>>var[366]>>var[367]>>var[368]>>var[369]>>var[370]>>var[371]>>var[372]>>var[373]>>var[374]>>var[375]>>var[376]>>var[377]>>var[378]>>var[379]>>var[380]>>var[381]>>var[382]>>var[383]>>var[384]>>var[385]>>var[386]>>var[387]>>var[388]>>var[389]>>var[390]>>var[391]>>var[392]>>var[393]>>var[394]>>var[395]>>var[396]>>var[397]>>var[398]>>var[399]>>var[400]>>
var[401]>>var[402]>>var[403]>>var[404]>>var[405]>>var[406]>>var[407]>>var[408]>>var[409]>>var[410]>>var[411]>>var[412]>>var[413]>>var[414]>>var[415]>>var[416]>>var[417]>>var[418]>>var[419]>>var[420]>>var[421]>>var[422]>>var[423]>>var[424]>>var[425]>>var[426]>>var[427]>>var[428]>>var[429]>>var[430]>>var[431]>>var[432]>>var[433]>>var[434]>>var[435]>>var[436]>>var[437]>>var[438]>>var[439]>>var[440]>>var[441]>>var[442]>>var[443]>>var[444]>>var[445]>>var[446]>>var[447]>>var[448]>>var[449]>>var[450]>>var[451]>>var[452]>>var[453]>>var[454]>>var[455]>>var[456]>>var[457]>>var[458]>>var[459]>>var[460]>>var[461]>>var[462]>>var[463]>>var[464]>>var[465]>>var[466]>>var[467]>>var[468]>>var[469]>>var[470]>>var[471]>>var[472]>>var[473]>>var[474]>>var[475]>>var[476]>>var[477]>>var[478]>>var[479]>>var[480]>>var[481]>>var[482]>>var[483]>>var[484]>>var[485]>>var[486]>>var[487]>>var[488]>>var[489]>>var[490]>>var[491]>>var[492]>>var[493]>>var[494]>>var[495]>>var[496]>>var[497]>>var[498]>>var[499]>>var[500]>>
var[501]>>var[502]>>var[503]>>var[504]>>var[505]>>var[506]>>var[507]>>var[508]>>var[509]>>var[510]>>var[511]>>var[512]>>var[513]>>var[514]>>var[515]>>var[516]>>var[517]>>var[518]>>var[519]>>var[520]>>var[521]>>var[522]>>var[523]>>var[524]>>var[525]>>var[526]>>var[527]>>var[528]>>var[529]>>var[530]>>var[531]>>var[532]>>var[533]>>var[534]>>var[535]>>var[536]>>var[537]>>var[538]>>var[539]>>var[540]>>var[541]>>var[442]>>var[543]>>var[544]>>var[545]>>var[546]>>var[547]>>var[548]>>var[549]>>var[550]>>var[551]>>var[552]>>var[553]>>var[554]>>var[555]>>var[556]>>var[557]>>var[558]>>var[559]>>var[560]>>var[561]>>var[562]>>var[563]>>var[564]>>var[565]>>var[566]>>var[567]>>var[568]>>var[569]>>var[570]>>var[571]>>var[572]>>var[573]>>var[574]>>var[575]>>var[576]>>var[577]>>var[578]>>var[579];

}
/*
   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] = 1;
   gMinuit->mnexcm("SET ERR", arglist ,4,ierflg);
   arglist[0] = 2;
   gMinuit->mnexcm("SET STR", arglist ,4,ierflg);
   
   //arglist[0] = 0.1;

//double par0 = gRandom->Uniform(0.65,0.75);
// double par1 = gRandom->Uniform(0.25,.35);
 //double par2 = gRandom->Uniform(-1.1,-1.001);
   //gMinuit->mnexcm("SET EPS", arglist, 1, ierflg); 
// Set starting values and step sizes for parameters
   gMinuit->mnparm(0, "omega_m", 0.3,0.01, 0,0,ierflg);
   gMinuit->mnparm(1, "wd", -1.1,0.0001, 0,0,ierflg);
   gMinuit->mnparm(2, "h   ", 0.7,0.01, 0,0,ierflg);
   //gMinuit->mnparm(3, "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("HESSE", 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(3,amin);
//printf(" \n chisq=%.10f",amin);
*/

}