//! Usage:
//! root
//! .L MinuitError.C++
//! MinuitError()

// Standard
#include <iostream>
#include <vector>
using namespace std;

// ROOT
#include "TH1.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TRandom.h"
#include "TROOT.h"
#include "TPluginManager.h"
#include "TStopwatch.h"

#include "Math/WrappedTF1.h"
#include "Math/WrappedMultiTF1.h"
#include "Fit/BinData.h"
#include "Fit/UnBinData.h"
#include "HFitInterface.h"
#include "Fit/Fitter.h"


/******************************************************************************/

class ParallelFitter
{
public:
  //! Default constructor
  ParallelFitter() { m_ThreadNextItem = 0; }
  //! Default destructor
  ~ParallelFitter() {};
  
  //! Analyze what eveer needs to be analyzed...
  bool Analyze();
 
  //! In this function the parallel calibration happens
  bool Calibrate();


private:
  //! Input: some histograms
  vector<TH1D*> m_Histograms;
  //! Output: Mean of the fitted Gaussian
  vector<double> m_Means;
  
  //! ID of the next item to be processed
  unsigned int m_ThreadNextItem;

};


/******************************************************************************
 * The fit function - Gaus
 */
class Gaus : public ROOT::Math::IParamFunction { 
public:
  void SetParameters(const double* p) { std::copy(p,p+NPar(),fp);}
  const double* Parameters() const { return fp; }
  ROOT::Math::IGenFunction* Clone() const { 
    Gaus* g = new Gaus(); 
    g->SetParameters(fp);
    return g;
  };
  unsigned int NPar() const { return 4; }
private:
  double DoEvalPar(double x, const double* par) const { 
    double fitval = par[0];
    double arg = 0;

    if (par[3] != 0) arg = (x - par[2])/par[3];
    fitval += par[1]*TMath::Exp(-0.5*arg*arg);

    return fitval;
  }
  double fp[4];  
};


/******************************************************************************
 * Do whatever analysis is necessary
 */
bool ParallelFitter::Analyze()
{
  // Sanity checks:
  // We will be using Minuit2 for fitting since it is thread safe, so check it is there
  if (gROOT->GetPluginManager()->FindHandler("ROOT::Math::Minimizer", "Minuit2") == 0)  {
    cout<<"We need Minuit2 for parallel fitting!"<<endl;
    return false;
  }
  
  
  // Create the list of histograms:
  unsigned int Histograms = 50;
  unsigned int Events = 5000;
  cout<<"Creating "<<Histograms<<" histograms with "<<Events<<" events for parallel fitting..."<<endl;
  for (unsigned int h = 0; h < Histograms; ++h) {
    TH1D* Hist = new TH1D("", "Hist", 100, -10, 10);
    for (unsigned int i = 0; i < Events; ++i) {
      double Val =  gRandom->Gaus(0, 1);
      Hist->Fill(Val);
    }
    m_Histograms.push_back(Hist);
    m_Means.push_back(-999);
  }
  TCanvas* HC = new TCanvas();
  HC->cd();
  m_Histograms[0]->Draw();
  HC->Update();
    
  Calibrate();
  
  // Create a histogram of all means
  TH1D* Means = new TH1D("", "Distribution of all fitted Gaussian means", 100, -0.1, 0.1);
  for (unsigned int i = 0; i < m_Means.size(); ++i) {
    //cout<<"Mean: "<<m_Means[i]<<endl;
    Means->Fill(m_Means[i]);  
  }
  
  TCanvas* C = new TCanvas();
  C->cd();
  Means->Draw();
  C->Update();
  C->SaveAs("Means.pdf");
  
  // Delete the histograms
  for (unsigned int i = 0; i < m_Histograms.size(); ++i) {
    delete m_Histograms[i]; 
  }
  m_Histograms.clear();
  m_Means.clear();
  
  return true;
}


/******************************************************************************
 * Perform the parallel calibration 
 */
bool ParallelFitter::Calibrate()
{
  Gaus FitFunction;
  ROOT::Fit::Fitter TheFitter; 
  TheFitter.Config().SetMinimizer("Minuit2");
  TheFitter.SetFunction(FitFunction);

  while (true) {
    unsigned int ID = m_ThreadNextItem;
    ++m_ThreadNextItem;
    
    if (ID >= m_Histograms.size()) break;
    
    ROOT::Fit::BinData d;
    d.Initialize(m_Histograms[ID]->GetNbinsX(), 1, ROOT::Fit::BinData::kValueError);
    for (int p = 1; p <= m_Histograms[ID]->GetNbinsX(); ++p) {
      d.Add(m_Histograms[ID]->GetBinCenter(p), m_Histograms[ID]->GetBinContent(p), sqrt(m_Histograms[ID]->GetBinContent(p)));
    }
    
    TheFitter.Config().ParSettings(0).SetName("Offset");
    TheFitter.Config().ParSettings(0).SetValue(0);
    TheFitter.Config().ParSettings(0).SetLimits(0, 0);
    TheFitter.Config().ParSettings(1).SetName("Height");
    TheFitter.Config().ParSettings(1).SetValue(1);
    TheFitter.Config().ParSettings(1).SetLimits(0, 10000);
    TheFitter.Config().ParSettings(2).SetName("Mean");
    TheFitter.Config().ParSettings(2).SetValue(1);
    TheFitter.Config().ParSettings(2).SetLimits(-10, 10);
    TheFitter.Config().ParSettings(3).SetName("Sigma");
    TheFitter.Config().ParSettings(3).SetValue(1);
    TheFitter.Config().ParSettings(3).SetLimits(0.1, 10);

    bool ReturnCode = false;
    ReturnCode = TheFitter.Fit(d);
   
    if (ReturnCode == true) {
      //cout<<"Fit: "<<TheFitter.Result().Parameter(0)<<":"<<TheFitter.Result().Parameter(1)<<":"<<TheFitter.Result().Parameter(2)<<":"<<TheFitter.Result().Parameter(3)<<endl;
      m_Means[ID] = TheFitter.Result().Parameter(2);
    }
  }
  
  return true;
}


/******************************************************************************
 * Main program
 */
int MinuitError()
{
  //int Threads = 4;
  
  ParallelFitter P;
  P.Analyze();

  cout<<"Program finished!"<<endl;

  return 0;
}