#include<iostream>
#if !defined(__CINT__) || defined(G__DICTIONARY)
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnUserFcn.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/ContoursError.h"
#endif
#include  "NeutronMultiFcn.h"
#include <vector>
#include <fstream>
#include <TCanvas.h>
#include <TStyle.h>
#include <TROOT.h>
#include <TAxis.h>
#include <TLine.h>
#include <TFile.h>
#include "TF1.h"
#include "TH1.h"

using namespace ROOT::Minuit2;

double nm(double *x, double* par)
{
  double xx = *x;
  NeutronMulti nm(par[0],par[1],par[2],par[3],par[4]);
  return nm(xx);
}


int DemoNMSim(){
  
  std::ifstream inpfile("F7.dat");
  // TH1D* h1 = new TH1D("h1","h2",100,0,20);
  std::vector<double> the_E;
  std::vector<double> the_Y;
  std::vector<double> the_Yer;
  int i=0;

  while(inpfile.good())
    { i++;
      double en =0;
      double yi =0;
      double yerr =0;
      inpfile >> en >> yi >> yerr ;
      the_E.push_back(en);
      the_Y.push_back(yi);
      the_Yer.push_back(yerr);
      cout<<i<<endl;
    }

  std::cout<<"***************"<<std::endl;
  NeutronMultiFcn nmfcn(the_E, the_Y, the_Yer);
  
  MnUserParameters upar;
  upar.Add("a1",2.0,0.2);
  upar.Add("a2",2.2,0.2);
  upar.Add("a3",2.3,0.2);
  upar.Add("t2",2,0.2);
  upar.Add("t3",3,0.3);

  //access parameters names to set limits
  // upar.SetLimits("a1",0.001,5);
  //upar.SetLimits("a2",0.001,5);
  //upar.SetLimits("a3",0.001,5);
  //upar.SetLimits("t2",0.001,5);
  //upar.SetLimits("t3",0.001,5);

  MnMigrad migrad(nmfcn,upar);
  migrad();
  FunctionMinimum  min = migrad();
  
  //  migrad.Fix("cons");

  MnUserParameterState fxn(upar);  
  double ta1=min.UserState().Value("a1"); 
  double ta2=min.UserState().Value("a2"); 
  double ta3=min.UserState().Value("a3");
  double tt2=min.UserState().Value("t2"); 
  double tt3=min.UserState().Value("t3");
  std::cout<<ta1<<"  "<<ta2<<std::endl;
  std::vector<double> par;
  par.push_back(ta1);par.push_back(ta2);par.push_back(ta3);
   par.push_back(tt2);par.push_back(tt3);

   cout<<min<<endl;
   // cout<<"The Chi2/ndf is: "<<nmfcn(par)<<"/"<<the_E.size()<<std::endl;
   // cout<<"The Chi2/ndf is: "<<nmfcn(par)/the_E.size()<<std::endl;

  // width=1.17;
  TF1* f1 = new TF1("ff",nm,2,20,10);
  f1->SetParameters(2.3,2.2,2.5,2,3);
  // myHist->GetYaxis()->SetRangeUser(1.,1e5);
    f1->Draw();
 
 
  return 0;  
}