// A stand-alone C++ program to illustrate the use of ROOT class
// TMinuit for function minimization.  

// Glen Cowan
// RHUL Physics
// 6 December 2004

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <string>
#include <vector>

#include <TMinuit.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TAxis.h>
#include <TLine.h>
#include <TPostScript.h>

using namespace std;

// Declare pointer to data as global (not very elegant)

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

// The pdf to be fitted, here an exponential. 
// First argument needs to be a pointer to plot with the TF1 class.

double expPdf(double* xPtr, double par[]){        
  double x = *xPtr;
  double xi = par[0];     // mean of x
  double f = 0;
  if ( xi > 0. ) { f = (1.0/xi) * exp(-x/xi); }
  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: ";
  getline(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();

}

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

void plotFit(TMinuit* minuit, TCanvas* canvas){

  int npar = minuit->GetNumPars();
  double outpar[npar], err[npar];
  for (int i=0; i<npar; i++){
    minuit->GetParameter(i,outpar[i],err[i]);
  }
  
  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);

  double xmin = 0.0;
  double xmax = 5.0;

  static TF1* func = 0;
  if ( func != 0 ) { delete func; }
  func = new TF1("funcplot", expPdf, xmin, xmax, npar);
  func->SetParameters(outpar);
  func->Draw();
  func->GetXaxis()->SetTitle("x");
  func->GetYaxis()->SetTitle("f(x;#xi)");

  vector<double> xVec = *xVecPtr;
  const double tickHeight = 0.1;
  TLine tick;
  for (int i=0; i<xVec.size(); i++){
    tick.DrawLine(xVec[i], 0, xVec[i], tickHeight);
  }
  canvas->Update();
  canvas->Show();

}

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

// fcn passes back f = -2*ln L; this is the function to be minimized.

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

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

// Compute log-likelihood function

  double lnL = 0.0;
  for (int i=0; i<n; i++){
    double x = xVec[i];
    double pdf = expPdf(&x, par);
    if ( pdf > 0.0 ) {
      lnL += log(pdf);    // need positive f
    }
  }
  f = -2.0 * lnL;         // factor of -2 so minuit gets the errors right
 
}                         // end of fcn

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

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

  TApplication* theApp = new TApplication("App", 0, 0);

// Initialize minuit, set initial values etc. of parameters.

  const int npar = 1;
  TMinuit* minuit = new TMinuit(npar);
  minuit->SetFCN(fcn);

  double par[npar]={1.00};
  double stepSize[npar]={0.1};
  double minVal[npar]={0.000001};
  double maxVal[npar]={100.00};
  string parName[npar]={"xi"};
  for (int i=0; i<npar; i++){
    minuit->DefineParameter(i, parName[i].c_str(), 
      par[i], stepSize[i], minVal[i], maxVal[i]);
  }

// Main user dialog

  bool keepGoing = true;
  TCanvas* canvas = new TCanvas();

  while ( keepGoing ) {

    cout << "Standard MINUIT command or in(put), plot, save, quit: ";
    cout << "Enter command: ";
    string command;
    getline(cin, command);

    if ( command == "quit" || command == "exit" ) {
      keepGoing = false;
    }
    else if ( command == "input" || command == "in" ) {
      getData(xVecPtr);
    }
    else if ( command == "plot" ) {
      plotFit(minuit, canvas);
    }
    else if ( command == "save" ) {
      cout << "output file name: ";
      string fileName;
      getline(cin, fileName);
      TPostScript file(fileName.c_str(), 113);     // 113 makes eps
      canvas->Draw();
      file.Close();
    }
    else {
      minuit->Command ( command.c_str() );
    }


  }         //     while keepGoing;

  delete minuit, xVecPtr, theApp;
  return 0;

}