/// \file
/// \ingroup tutorial_fit
/// \notebook
/// 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 Lorenzo Moneta

#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
using namespace std;


const int N = 2; // number of fit functions
// definition of shared parameter
// background function
int iparB[2] = { 0,      // exp amplitude in B histo
                 2    // exp common parameter
};

// signal + background function
int iparSB[5] = { 1, // exp amplitude in S+B histo
                  2, // exp common parameter
                  3, // gaussian amplitude
                  4, // gaussian mean
                  5  // gaussian sigma
};

// Create the GlobalChi2 structure

struct GlobalChi2 {
  GlobalChi2(  std::vector<ROOT::Fit::Chi2Function> &f)
   {
     for(int k=0; k<N; k++) fChi2.push_back(&f[k]);
   }
   // parameter vector is first background (in common 1 and 2)
   // and then is signal (only in 2)
  double operator() (const double *par) const
   {
     // use std::vector in case N is not known at compiled time
      std::vector<std::vector<double>> p(N);
        p[0].resize(2); p[1].resize(5);
        for (int i = 0; i < 2; ++i) p[0][i] = par[iparB[i] ];
        for (int i = 0; i < 5; ++i) p[1][i] = par[iparSB[i] ];
      double sum = 0; 
      for(int k=0; k<N; k++) sum += (*fChi2[k])(p[k].data());
      return sum;  
   }
  std::vector< ROOT::Math::IMultiGenFunction *> fChi2;
};

void combinedFitMod3() { //
   TH1D *h[N];
   h[0] = new TH1D("hB","histo B",100,0,100);
   h[1] = new TH1D("hSB","histo S+B",100, 0,100);
   TF1 *f[N];
   f[0] = new TF1("fB","expo",0,100);
   f[0]->SetParameters(1,-0.05);
   TF1 *fS = new TF1("fS","gaus",0,100);
   fS->SetParameters(1,30,5);

   h[0]->FillRandom("fB");
   h[1]->FillRandom("fB",2000);
   h[1]->FillRandom("fS",1000);

   // perform now global fit

   f[1] = new TF1("fSB","expo + gaus(2)",0,100);

   ROOT::Fit::DataOptions opt;
   std::vector<ROOT::Fit::DataRange> range;
   ROOT::Fit::BinData data[N];             // working
   //std::vector<ROOT::Fit::BinData> data; //not working
   std::vector< ROOT::Math::WrappedMultiTF1 > wf;
   std::vector<ROOT::Fit::Chi2Function> chi2;
   for (int i = 0; i < N; ++i) {
     range.push_back(ROOT::Fit::DataRange());
     range[i].SetRange(10,90-i*40);
     data[i] = ROOT::Fit::BinData(opt, range[i]);
     ROOT::Fit::FillData(data[i], h[i]);                 // working
     //data.push_back(ROOT::Fit::BinData(opt,range[i])); //not working
     wf.push_back(ROOT::Math::WrappedMultiTF1(*f[i],1));
     chi2.push_back(ROOT::Fit::Chi2Function(data[i], wf[i]));
   }
   GlobalChi2 globalChi2(chi2);

   ROOT::Fit::Fitter fitter;

   const int Npar = 6;
   double par0[Npar] = { 5,5,-0.1,100, 30,10};

   // create before the parameter settings in order to fix or set range on them
   fitter.Config().SetParamsSettings(6,par0);
   // fix 5-th parameter
   fitter.Config().ParSettings(4).Fix();
   // set limits on the third and 4-th parameter
   fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4);
   fitter.Config().ParSettings(3).SetLimits(0,10000);
   fitter.Config().ParSettings(3).SetStepSize(5);

   fitter.Config().MinimizerOptions().SetPrintLevel(0);
   fitter.Config().SetMinimizer("Minuit2","Migrad");

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

   unsigned int dataSize = 0;
   for (int i=0; i<N; i++) dataSize += data[i].Size();
   fitter.FitFCN(6,globalChi2,0,dataSize,true);
   
   ROOT::Fit::FitResult result = fitter.Result();
   result.Print(std::cout);

   TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms",10,10,700,700);
   c1->Divide(1,N);
   gStyle->SetOptFit(1111);
   for (int i = 0; i < N; ++i) {
     c1->cd(i+1);
     f[i]->SetFitResult( result, iparB);
     if(i>0) f[i]->SetFitResult( result, iparSB);
     f[i]->SetLineColor(4-2*i);
     h[i]->Draw();
     f[i]->Draw("same");
   }

}