// Working code adapted from combinedFit_LL.C (by L. Moneta - Dec 2010)
// 
// ---------------------------------------------------------------------------------------
//+ Combined (simultaneous) fit of two histogram with separate functions 
//  and some common parameters
//
// See http://root.cern.ch/phpBB3//viewtopic.php?f=3&t=11740#p50908
// for a modified version working with Fumili or GSLMultiFit 
//
// N.B. this macro must be compiled with ACliC 
//
//Author: L. Moneta - Dec 2010
// ---------------------------------------------------------------------------------------


#include "Fit/Fitter.h"
#include "Fit/BinData.h"
//#include "Fit/Chi2FCN.h"
#include "Fit/PoissonLikelihoodFCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"



void propagateSystematicErrors() { 

  TH1D * h = new TH1D("h","Exponential decay",100,0,1000);

  TF1 * f = new TF1("f","expo",h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
  f->SetParameters(5,-0.01);
  h->FillRandom("f");

  ROOT::Math::WrappedMultiTF1 wf(*f,1);

  ROOT::Fit::DataOptions opt; 
    opt.fUseEmpty = true;
    opt.fExpErrors = true;
  ROOT::Fit::DataRange rangeB; 
  // set the data range
  rangeB.SetRange( h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
  
  ROOT::Fit::BinData data_h(opt,rangeB); 
  
  ROOT::Fit::FillData(data_h, h);

  ROOT::Fit::PoissonLLFunction likelihood_F(data_h, wf);

  ROOT::Fit::Fitter fitter;

  // create before the parameter settings in order to fix or set range on them
  fitter.Config().SetParamsSettings( f->GetNpar(), f->GetParameters() );
  // fix the tau-decay factor
  //fitter.Config().ParSettings(1).Fix();
  
  fitter.Config().MinimizerOptions().SetPrintLevel(0);
  fitter.Config().SetMinimizer("Minuit2","Migrad"); 

  // fit FCN function directly 
  // (specify optionally data size and flag to indicate if it is a chi2 fit)

  double meanPar = -0.01;
  double sigmaPar = 0.001; 
  auto nll = [&](const double *x) { return likelihood_F(x) - std::log( ROOT::Math::normal_pdf(x[1],sigmaPar, meanPar) ); };
  
  fitter.FitFCN(2,nll,0,data_h.Size(),false);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  // Draw fit
  TCanvas * c1 = new TCanvas("c1","",10,10,1280,800);
    gStyle->SetOptFit(1111);
    f->SetFitResult( result );
    f->SetLineColor(kBlue);
    h->GetListOfFunctions()->Add(f);
    h->Draw(); 
  
  return;
  
}