#define  myntupme_cxx
#include "myntupme.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TLorentzVector.h>
#include <iostream>
#include <TGraphErrors.h>
#include <vector>
#include "TVector3.h"
#include "TPostScript.h"
#include "TF1.h"
#include "TRandom.h"
#include "TCutG.h"
#include <TLatex.h>

#include "TMath.h"
#include "TLorentzVector.h"
#include "Riostream.h"
#include "TFitter.h"
#include "TH1.h"
#include "TH2.h"
#include "TMinuit.h"
#include "TGenPhaseSpace.h"
#include "TRandom2.h"
#include "TROOT.h"
#include "TFile.h"

#include <cstdlib>                                      // C standard library
#include <cstdio>                                       // C I/O (for sscanf)
#include <cstring>                                      // string manipulation
#include <fstream>                                      // file I/O



double a,b,c,d,gout,h, par[5];
int nDim, flg;
Bool_t         debug    = false; // Debug flag for some printouts.

void chisqFunction(int& nDim, double* gout, double& a, double par[], int flg);
bool KinematicFit(TFitter *minimizer, double *th, double *ph, double *e);

double a=0;

void chisqFunction(int& nDim, double* gout, double& a , double par[], int flg)
{
 h = ((par[0]-thetaProtons[0])/phiProtons[0])*((par[0]-thetaProtons[0])/phiProtons[0]); 
a=(a+h);
 
}


bool KinematicFit(TFitter *minimizer, double *th, double *ph, double *e)
{
  bool converged = false;
  int count = 0;
  Int_t ierr;
 
  double line_points[4];
  double result_point[2];

  thetaProtons[0] = th[0];
  phiProtons[0] = ph[0];

 cout << "ok1" <<"\n";

  while (!converged)
    {
cout << "ok2" <<"\n";	
      // First solve for proton 2 fixed
      
      minimizer->SetParameter(0,       // variable number
			      "e1",    // variable name
			      e[0],    // variable value
			      phiProtons[0], // range to search for solution
			      0,
			      0);
cout << "ok3" <<"\n";

      // Keep E1 as free parameter and fix E2 
      minimizer->ReleaseParameter(0);
//      minimizer->FixParameter(1);      
cout << "ok4" <<"\n";	
      ierr = minimizer->ExecuteCommand("MIGRAD",0,0); // Find solution for E1
cout << "ok5" <<"\n";
      if (ierr) return false;
      
      line_points[0] = minimizer->GetParameter(0);

	  cout << "1-Fitter E1: " << count << "/" << e[0] << "-->" << line_points[0] << endl;

}
}	

void myntupme::Loop()
{

double chi,chisq,chisqnew,dsnew,ds3ncnew;
double chisqo,chisqnewo;
chisqo=0;
chisqnewo=0;
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;

   double chi,ph,chi1,chi2,chix,chix1,chix2,chit,chit1,chit2,chit3,chi3,chitx,chitx1,chitx2,chitx3,chix3,t1[20],t2[20],chics,chics1,chics2,chicst,chicst1,chicst2,chicst3,chics3;

   double chis,chis1,chis2,chixs,chixs1,chixs2,chis3,chixs3,chicss,chicss1,chicss2,chicss3;
   double chisp,chisp1,chisp2,chixsp,chixsp1,chixsp2,chisp3,chixsp3,chicssp,chicssp1,chicssp2,chicssp3;
double chism,chism1,chism2,chixsm,chixsm1,chixsm2,chism3,chixsm3,chicssm,chicssm1,chicssm2,chicssm3;

   TGraphErrors *gcs2=new TGraphErrors();
   TGraphErrors *gcs1=new TGraphErrors();

int sm,h,k,q,cs,cst;


ifstream Sinput;

Sinput.open("conf.dat");
  if(Sinput) fprintf(stdout,"file  conf is open.\n ");
  else {
    fprintf(stdout,"can't open file conf \n");
  }

for(k=0; k<10; k++) { 
  
    Sinput >> t1[k]  >> t2[k];
}
Sinput.close();

for (int m=0; m<1; m++){ 
cs=0;
sm=0;
ph=m*20+20;
chi=0, chi2=0, chi3=0, chix=0, chix1=0, chix2=0, chix3=0, chics=0, chics1=0, chics2=0 ,chics3=0;
chis2=0,chixs2=0,chis3=0,chixs3=0,chicss2=0,chicss3=0;
chisp=0,chisp1=0,chisp2=0,chixsp=0,chixsp1=0,chixsp2=0,chisp3=0,chixsp3=0,chicssp=0,chicssp1=0,chicssp2=0,chicssp3=0; chism=0,chism1=0,chism2=0,chixsm=0,chixsm1=0,chixsm2=0,chism3=0,chixsm3=0,chicssm=0,chicssm1=0,chicssm2=0,chicssm3=0;



  bool converged = false;
  int count = 0;
  int countt = 0;
  Int_t ierr;
 

  TFitter* minimizer = new TFitter(2);
  minimizer->GetMinuit()->SetPrintLevel(-1); // Be silent!
  
  // Tell the minimizer about the function to be minimzed
  minimizer->SetFCN(chisqFunction);

 
   for (Long64_t jentry=0; jentry<2;jentry++) {
     Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

ds=ds*1000;
esysds=esysds*1000;
ds3nc=ds3nc*1000;
ds3n=ds3n*1000;
eds=eds*1000;
dsnn=dsnn*1000;
ds_tm=ds_tm*1000;

/*    energy_protons = ds;
    theta_protons = ds3nc;
    phi_protons = eds;
*/
   bool success = KinematicFit(minimizer,ds3nc, eds, ds);

}
}
}