#define znevzat_cxx
#include "TROOT.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TStopwatch.h"
#include "TTree.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TChain.h"
#include "TStyle.h"
#include "TMath.h"
#include "TFormula.h"
#include "TProfile.h"
#include "TF1.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
//gSystem->Load("libMinuit");
#include "TMinuit.h"

#include "znevzat.h"


//using namespace std;

//int main(){
//  znevzat r;
//  r.Loop();
//}
const Int_t ndata=20000;
Double_t deltaPHI[ndata];
Double_t deltaZ[ndata];
Double_t errorDZ[ndata];
Double_t deltaTHETA[ndata];
Double_t sigmaDZ[300];
Double_t thetabin[50];
Double_t mombin[50];
Double_t rad2deg=180.0/3.14159;

//______________________________________________________________________________
Double_t func(Double_t Df,Double_t *par)
{
 Double_t value = (par[0]*Df+par[1]);
 return value;
}

//============================================================================================
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag, Int_t ndata)
{
	Int_t i;

////calculate chisquare
	
	Double_t chisq=0;
	Double_t chi;
	for(i=0;i<ndata;i++) {
		chi=(deltaZ[i]-func(deltaPHI[i],par));
		chisq+=chi*chi;
	}	
   f = chisq;
}
//=============================================================================================


void znevzat::Loop()
{
   if (fChain == 0) return;
   Int_t nentries = Int_t(fChain->GetEntriesFast());

//====================== DEFINE THINGS===============================================
TH2D *hh = new TH2D("hh", "temporary histogram", 100, -1., 1.,100,-1.,1.);

for(Int_t j=0;j<ndata;j++){
deltaPHI[j]=deltaZ[j]=errorDZ[j]=deltaTHETA[j]=0.0;}
//==============================================================================

   Int_t nbytes = 0, nb = 0;
   //for (Int_t jentry=0; jentry<nentries;jentry++) {
   for (Int_t jentry=0; jentry<ndata;jentry++) {
//====================== FILL ARRAYS and HISTOGRAMS =============================

	if(npro==1&&abs(prz-gprz)<1.0&&abs(prf-gprf)<0.5&&abs(gprm-0.95)<0.1&&abs(57.*prt-35.)<10.){
	hh->Fill((prz-gprz)/prm,prf-gprf);
//	hh->Fill(gprt*rad2deg,gprm);
	}
	
	deltaPHI[jentry]=prf-gprf;
	deltaZ[jentry]=prz-gprz;
	deltaTHETA[jentry]=prt-gprt;
	errorDZ[jentry]=prz-rprz;
	
//===============================================================================
      Int_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
   }
//====================== ANALYZE RESULTS ========================================


TCanvas *c1 = new TCanvas("c1","",0,0,600,400);
gStyle->SetOptFit(0111);
gPad->SetGrid();

hh->Draw();
TProfile *prof = hh->ProfileX();
prof->SetLineColor(kRed);
prof->Draw("same");
prof->Draw("");
prof->Fit("pol1","R","",-0.6,0.6); //<-----

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

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

// Set starting values and step sizes for parameters
   static Double_t vstart[2] = {0.0 , 0.0};
   static Double_t step[2] = {0.0001 , 0.0001};
   gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
   gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);

// Now ready for minimization step
   arglist[0]=100000000;
   arglist[1]=0.001;  //<-----
   gMinuit->mnexcm("MIGRAD", 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(3,amin);
}