#include "Minuit2/GaussDataGen.h"
#include "Minuit2/FCNBase.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MinimumPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/VariableMetricMinimizer.h"
#include "Minuit2/MnUserParameters.h"

#include "TSystem.h"

#include <iostream>


class GaussFunction {

public:

  GaussFunction(double mean, double sig, double constant) :
    theMean(mean), theSigma(sig), theConstant(constant) {}

  ~GaussFunction() {}

  double m() const {return theMean;}
  double s() const {return theSigma;}
  double c() const {return theConstant;}

  double operator()(double x) const {

    return
       c()*exp(-0.5*(x-m())*(x-m())/(s()*s()))/(sqrt(2.*M_PI)*s());
  }

private:

  double theMean;
  double theSigma;
  double theConstant;
};






class GaussFcn : public ROOT::Minuit2::FCNBase {

public:

  GaussFcn(const std::vector<double>& meas,
           const std::vector<double>& pos,
     const std::vector<double>& mvar) : theMeasurements(meas),
                      thePositions(pos),
                theMVariances(mvar),
                theErrorDef(1.) {}

  ~GaussFcn() {}

  double Up() const {return theErrorDef;}
  //virtual double operator()(const std::vector<double>&) const;
  
  double operator()(const std::vector<double>& par) const {

  assert(par.size() == 3);
  GaussFunction gauss(par[0], par[1], par[2]);

  double chi2 = 0.;
  for(unsigned int n = 0; n < theMeasurements.size(); n++) {
    chi2 += ((gauss(thePositions[n]) - theMeasurements[n]) *
             (gauss(thePositions[n]) - theMeasurements[n]) /
             theMVariances[n]);
  }

  return chi2;
  }

  std::vector<double> measurements() const {return theMeasurements;}
  std::vector<double> positions() const {return thePositions;}
  std::vector<double> variances() const {return theMVariances;}

  void setErrorDef(double def) {theErrorDef = def;}

private:


  std::vector<double> theMeasurements;
  std::vector<double> thePositions;
  std::vector<double> theMVariances;
  double theErrorDef;
};







using namespace std;

int TestMinuit() {
	
	gSystem->Load("libMinuit2.so");
	
  // generate the data (100 data points)
  ROOT::Minuit2::GaussDataGen gdg(100);

  std::vector<double> pos = gdg.positions();
  std::vector<double> meas = gdg.measurements();
  std::vector<double> var = gdg.variances();


  // create FCN function
  GaussFcn theFCN(meas, pos, var);

  // create initial starting values for parameters
  double x = 0.;
  double x2 = 0.;
  double norm = 0.;
  double dx = pos[1]-pos[0];
  double area = 0.;
  for(unsigned int i = 0; i < meas.size(); i++) {
    norm += meas[i];
    x += (meas[i]*pos[i]);
    x2 += (meas[i]*pos[i]*pos[i]);
    area += dx*meas[i];
  }
  double mean = x/norm;
  double rms2 = x2/norm - mean*mean;
  double rms = rms2 > 0. ? sqrt(rms2) : 1.;

  {
    // demonstrate minimal required interface for minimization
    // create Minuit parameters without names

    // starting values for parameters
    std::vector<double> init_par;
    init_par.push_back(mean);
    init_par.push_back(rms);
    init_par.push_back(area);

    // starting values for initial uncertainties
    std::vector<double> init_err;
    init_err.push_back(0.1);
    init_err.push_back(0.1);
    init_err.push_back(0.1);

    // create minimizer (default constructor)
    ROOT::Minuit2::VariableMetricMinimizer theMinimizer;

    // minimize
    ROOT::Minuit2::FunctionMinimum min = theMinimizer.Minimize(theFCN, init_par, init_err);

    // output
    std::cout<<"minimum: "<<min<<std::endl;
    }

	{
    // demonstrate standard minimization using MIGRAD
    // create Minuit parameters with names
    ROOT::Minuit2::MnUserParameters upar;
    upar.Add("mean", mean, 0.1);
    upar.Add("sigma", rms, 0.1);
    upar.Add("area", area, 0.1);

    // create MIGRAD minimizer
    ROOT::Minuit2::MnMigrad migrad(theFCN, upar);

    // minimize
    ROOT::Minuit2::FunctionMinimum min = migrad();

    // output
    std::cout<<"minimum: "<<min<<std::endl;
    }
	

  return 0;
};