#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <cmath>

#include "PLFit.h"
#include "PLUtil.h"
#include "fitPar.h"
#include <TRandom3.h>
#include <TFile.h>
#include <TNtuple.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TH1D.h>
#include <TMinuit.h>
#include <TGraph.h>
#include <TPostScript.h>

PLFit* plfGlobal;  
bool debug;

int fitPar (PLFit* plf, int npar, double& chi2Min, vector<double>& theta,
            vector<vector<double> >& V){

  plfGlobal = plf;                  // communicate to fcn via global

  const int maxpar = 100;
  double par[maxpar];
  string parName[maxpar];
  double stepSize[maxpar];
  double minVal[maxpar];
  double maxVal[maxpar];

  for (int i=0; i<npar; i++){
    par[i] = theta[i];
  }

  int status = 0;
  if ( npar == 0 ){          // just compute chi2, no fit

    double* deriv;
    int flag;
    fcn(npar, deriv, chi2Min, par, flag);

  }
  else {
 
// set up start values, step sizes, etc. for fit

    for (int i=0; i<npar; i++){    // some of these overridden below
      stepSize[i] = 0.1;
      minVal[i] = 0.;
      maxVal[i] = 0.;
    }  

    for (int i=0; i<npar; i++){
      stringstream ss;
      string parNumber;
      ss << i;
      ss >> parNumber;
      parName[i] = "p" + parNumber;
    }

    for (int i=0; i<npar; i++){
      if ( par[i] != 0 ) {
        stepSize[i] = abs( par[i] ) * 0.1;
      }
      else {
        stepSize[i] = 0.1;
      }
    }

// Create the TMinuit object and define parameters

    TMinuit* minuit = new TMinuit(npar);

    minuit->SetPrintLevel(-1);
    int ierr = 0;
    minuit->mnexcm("SET NOWarnings",0,0,ierr);

    minuit->SetFCN(fcn);
    for (int i=0; i<npar; i++){
      minuit->DefineParameter(i, parName[i].c_str(), 
      par[i], stepSize[i], minVal[i], maxVal[i]);
    }

// do the fit and extract results

    status = minuit->Command("MIGRAD");
    if ( status != 0 ) {
      debug = true;
      status = minuit->Command("SIMPLEX");
      status = minuit->Command("MIGRAD");
      debug = false;
    }

    theta.clear();
    int nFreePar = minuit->GetNumFreePars();
    double fitpar[nFreePar], err[nFreePar];
    for (int i=0; i<npar; i++){
      minuit->GetParameter(i, fitpar[i], err[i]);
      theta.push_back(fitpar[i]);
    }

    double covMatrix[npar][npar];
    minuit->mnemat(&covMatrix[0][0], npar);
    V.clear();
    for (int i=0; i<npar; i++){
      vector<double> rowI;
      for (int j=0; j<npar; j++){
        rowI.push_back(covMatrix[i][j]);
      }
      V.push_back(rowI);
    }

    double fedm, errdef;
    int npari, nparx, istat;
    minuit->mnstat(chi2Min, fedm, errdef, npari, nparx, istat);

    delete minuit;

  }                  // only chi2 or do fit

  return status;     // 0 = normal completion, 4 = MIGRAD not converged

}

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

// fcn must be non-member function, uses global PLFit object fit

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

  PLFit* plf = plfGlobal;

  vector<double> pdfPar;            // for now not used
  vector<double> modelBeta;
  for (int i=0; i<npar; i++){
    modelBeta.push_back(par[i]);
  }
  double nuTot = plf->nuTot();
  int opt = plf->fitOption();
  plf->setModelHist(pdfPar, modelBeta);

  double chi2 = 0;
  for (int i=1; i<=plf->nBin(); i++){ 
    double n = plf->dataHist()->GetBinContent(i);
    double nu = plf->modelHist()->GetBinContent(i);
    if ( n > 0. && nu > 0. ) {
      chi2 += 2 * (n * log(n/nu) + nu - n);
    }
    else {
      chi2 += nu - n;
    }
  }
 
  f = chi2;     // pass back by reference

}