#include <stdlib.h>
#include <stdio.h>
#include <Riostream.h>
#include <vector>
#include <list>
#include <string>

#include <TSystem.h>
#include <TStyle.h>
#include <TROOT.h>

#include <TCanvas.h>
#include <TF1.h>
#include <TList.h>
#include <TString.h>
#include <TObject.h>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TMinuit.h>
#include <TRandom3.h>
#include <TGraph.h>
#include <TMath.h>
#include <TArc.h>


   TGraph *IO = new TGraph("Experimental_Elastic.txt");

   void fcn(int &npar, double *gin, double &f, double *x, int iflag)
   {
      npar = 6;

      //this creates a Test.txt file which is the input to the fortran executable in the next line. 
      gROOT->ProcessLine(".L creating_txt_for_chuck.C+");
      gROOT->ProcessLine(Form("creating_txt_for_chuck(%0.1f,%0.3f, %0.3f, %0.1f, %0.3f, %0.3f)",x[0],x[1],x[2],x[3],x[4],x[5]));
      
      //system call to run the fortran executable
      system("./chuck <Test.txt> Elastic.txt"); //Elastic.txt is the output of the Fortran code
 
      TGraph *gr_chuck = new TGraph("Elastic.txt"); //creating a graph with the output of the Fortran code

     

      //calculate residual 
      Double_t residual = 0;
      Double_t delta;
      //for (i=0;i<nbins; i++) {
      for(double i=0.0;i<10.0;i=i+0.2){
           delta  = (gr_chuck->Eval(i)-IO->Eval(i));
           residual += delta*delta;
      }
           f = residual;
    }

    

    void optical_potential_parameter(){

    Double_t value1,err1;
    Double_t value2,err2;
    Double_t value3,err3;
    Double_t value4,err4;
    Double_t value5,err5;
    Double_t value6,err6;
    
    Double_t vstart[6] = {-230.0, 1.495, 0.418, -150.0, 1.212, 0.418};

    TMinuit *gMinuit = new TMinuit(6);  //initialize TMinuit with a maximum of 6 params   

    gMinuit->SetFCN(fcn);
    
    gMinuit->DefineParameter(0, "VR", vstart[0], 5.0,  -400.0, 0.0);
    gMinuit->DefineParameter(1, "RR", vstart[1], 2.0,   0.0, 10.0);
    gMinuit->DefineParameter(2, "AR", vstart[2], 0.5, 0.0, 1.0);
    gMinuit->DefineParameter(3, "VI", vstart[3], 5.0,  -300.0, 0.0);
    gMinuit->DefineParameter(4, "RI", vstart[4], 2.0,   0.0, 10.0);
    gMinuit->DefineParameter(5, "AI", vstart[5], 0.5, 0.0, 1.0);

    //gMinuit->DefineParameter(0, "VR", vstart[0], 50,  0, 0);
    //gMinuit->DefineParameter(1, "RR", vstart[1], 2,   0, 0);
    //gMinuit->DefineParameter(2, "AR", vstart[2], 0.5, 0, 0);
    //gMinuit->DefineParameter(3, "VI", vstart[3], 50,  0, 0);
    //gMinuit->DefineParameter(4, "RI", vstart[4], 2,   0, 0);
    //gMinuit->DefineParameter(5, "AI", vstart[5], 0.5, 0, 0);

    gMinuit->Command("MIGRAD");
    gMinuit->Command("HESSE");
    gMinuit->Command("MINOS");
    gMinuit->Command("CALL FCN");     

    gMinuit->GetParameter(0, value1, err1);
    gMinuit->GetParameter(1, value2, err2);
    gMinuit->GetParameter(2, value3, err3);
    gMinuit->GetParameter(3, value4, err4);
    gMinuit->GetParameter(4, value5, err5);
    gMinuit->GetParameter(5, value6, err6);

    cout << "First parameter" << value1 << endl;


}