#include "toolKineFit.h"
#include "TROOT.h"
#include "TString.h"

// *****
toolKineFit::toolKineFit() {
  // set global constants
  gl_n = 0;
  gl_nMax = 0;
  // maximum number of parameters
  int maxpar = 6;
  // initialize minuit, which is global
  gMinuit = new TMinuit(maxpar);
  gMinuit->SetFCN(gl_fcn);  // gl_fcn is global
}

// *****
toolKineFit::~toolKineFit() {}

// *****
void
toolKineFit::setJets(std::vector<TLorentzVector>& jets) {
  // gl_jets is global
  gl_jets->clear();
  for (int j = 0; j < (int)jets.size(); ++j) gl_jets->push_back(jets[j]);
}

// *****
void
toolKineFit::setPrintLevel(int level, int ncalls) {
  gl_nMax = ncalls;
  int ierflg = 0;
  double arglist[1];
  arglist[0] = level;
  gMinuit->mnexcm("SET PRI", arglist, 1, ierflg);
}

// *****
bool
toolKineFit::fit() {
  // define starting values and setp sizes for parameters

  // get the data
  std::vector<TLorentzVector> jets = *gl_jets;  // gl_jets is global

  // initialize parameters
  // energies

  // get the error on the jet energy
  std::vector<double> eSigmas;
  gl_errors(jets, eSigmas);

  for (int i = 0; i < (int)jets.size(); ++i) {
    int ierflg = 0;
    // guess initial value
    double start = jets[i].E();
    // use guess initial error on the jet energy
    double step = eSigmas[i];
    // if minVal = maxVal = 0 then the parameter is unbounded
    double minVal = 0.;
    double maxVal = 1000000.;
    TString t = "e"; t += i;
    ierflg = 0;
    gMinuit->mnparm(i, t, start, step, minVal, maxVal, ierflg);
    // std::cout << "mnparm " << i << " ierflg = " << ierflg << std::endl;
  }

  // chi2 "up" parameter for error definition
  gMinuit->SetErrorDef(1);
  // start minimization
  gMinuit->Migrad();
  //
  return true;
}

// *****
bool
toolKineFit::isValid() {
  double fmin, fedm, errdef;
  int nvpar, nparx, icstat;
  gMinuit->mnstat(fmin, fedm, errdef, nvpar, nparx, icstat);
  // std::cout << "icstat = " << icstat << std::endl;
  return icstat==3;
}
  

// *****
std::vector<TLorentzVector>
toolKineFit::getJetsFit() {
  // call function with termination flag to fill vector of jetFits
  int ierflg = 0;
  double arglist[1] = {3};
  gMinuit->mnexcm("CAL1", arglist, 1, ierflg);
  return *gl_jetsFit;
}

// *****
double
toolKineFit::getFcnMin() {
  double fmin, fedm, errdef;
  int nvpar, nparx, icstat;
  gMinuit->mnstat(fmin, fedm, errdef, nvpar, nparx, icstat);
  return fmin;
}