//Test to find parameter a and b of a linear function y=ax+b
//g++ `root-config --cflags --glibs` -lMinuit affine.cpp 
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <vector>
#include <iomanip>
#include <time.h>



#include "/usr/local/root/include/TH1F.h"
#include "/usr/local/root/include/TFile.h"
#include "/usr/local/root/include/cfortran.h"
#include "/usr/local/root/include/TMinuit.h"
using namespace std;

vector<double>* xVecPtr = new vector<double>();
vector<double>* yVecPtr = new vector<double>();

//my function to minimize

double affine(double* xPtr, double par[]){        
  double x = *xPtr;
  double f = 0;
  f = par[0]*x + par[1];
  return f;
}

//-------------------------------------------------------------------------

//  function to read in the data from a file

void getData(vector<double>* xVecPtr){

  string infile;
  cout << "Enter name of input data file: ";
  cin >> infile;

  ifstream f;
  f.open(infile.c_str());
  if ( f.fail() ){
    cout << "Sorry, couldn't open file !" << endl;
    exit(1);
  }

  double x ;
  bool acceptInput = true;
  while ( acceptInput ) {
    f >> x;
    acceptInput = !f.eof();
    if ( acceptInput ) {
      xVecPtr->push_back(x);
    }
  }
  f.close();
  return;
}

void getData1(vector<double>* yVecPtr){

  string infile;
  cout << "Enter name of input data file: ";
  cin >> infile;

  ifstream f;
  f.open(infile.c_str());
  if ( f.fail() ){
    cout << "Sorry, couldn't open file !" << endl;
    exit(1);
  }

  double y ;
  bool acceptInput = true;
  while ( acceptInput ) {
    f >> y;
    acceptInput = !f.eof();
    if ( acceptInput ) {
      yVecPtr->push_back(y);
    }
  }
  f.close();
  return;
}


//-------------------------------------------------------------------------


void fcn(int& npar, double* deriv, double& f, double par[], int flag){

  vector<double> xVec = *xVecPtr; 
  int n = xVec.size();

  double lnL = 0.0;
  for (int i=0; i<n; i++){
    double x = xVec[i];
    double pdf = affine(&x, par);
    f = pdf;
  }
  return;
}

//-------------------------------------------------------------------------

int main(int argc, char **argv) {

  getData(xVecPtr);   
  getData1(yVecPtr);   

// Initialize minuit, set initial values etc. of parameters.
  TMinuit *gMinuit = new TMinuit(2);  
  gMinuit->SetFCN(fcn);

  double arglist[10];
  int error_flag = 0;

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

  gMinuit->mninit(5,6,7); 
  double f_null=0.0;
  error_flag = 0;

  gMinuit->mnparm(0,   "a",0.,0.1,f_null,f_null,error_flag);
  gMinuit->mnparm(1,  "b", 0., 0.1,f_null,f_null,error_flag);

  // minimize
  gMinuit->mnexcm("SIMPLEX",0,0,error_flag);
  // get parameters
  double pval[2],perr[2],plo[2],phi[2];
  TString para0, para1;
  int istat = 0;
  gMinuit->mnpout(0,para0,pval[0],perr[0],plo[0],phi[0],istat);
  gMinuit->mnpout(1,para1,pval[1],perr[1],plo[1],phi[1],istat);

  // display results
  cout << setw(8) << pval[0] << setw(8) << perr[0] << endl;

  return 0;

}