// combined fit of two histogram with separate functions

#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"


struct GlobalChi2 { 
   GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,  ROOT::Math::IMultiGenFunction & f2) : 
      fChi2_1(f1), fChi2_2(f2) {}

   // parameter vector is first background (in common 1 and 2) and then is signal (only in 2)
   double operator() (const double *par) const {
      double p1[2]; 
      p1[0] = par[0]; // exp amplitude in B histo
      p1[1] = par[2]; // exp common parameter 

      double p2[5]; 
      p2[0] = par[1]; // exp amplitude in S+B histo
      p2[1] = par[2]; // exp common parameter
      p2[2] = par[3]; // gaussian amplitude
      p2[3] = par[4]; // gaussian mean
      p2[4] = par[5]; // gaussian sigma
      return fChi2_1(p1) + fChi2_2(p2);
   } 

   const  ROOT::Math::IMultiGenFunction & fChi2_1;
   const  ROOT::Math::IMultiGenFunction & fChi2_2;
};

void combinedFit() { 

  TH1D * hB = new TH1D("hB","histo B",100,0,100);
  TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);

  TF1 * fB = new TF1("fB","expo",0,100);
  fB->SetParameters(1,-0.05);
  hB->FillRandom("fB");

  TF1 * fS = new TF1("fS","gaus",0,100);
  fS->SetParameters(1,30,5);

  hSB->FillRandom("fB",2000);
  hSB->FillRandom("fS",1000);

  // perform now global fit

  TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);

  ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
  ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);

  ROOT::Fit::BinData dataB; 
  ROOT::Fit::FillData(dataB, hB);

  ROOT::Fit::BinData dataSB; 
  ROOT::Fit::FillData(dataSB, hSB);

  ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
  ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);

  GlobalChi2 globalChi2(chi2_B, chi2_SB);

  ROOT::Fit::Fitter fitter;

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

  fitter.FitFCN(6,globalChi2,par0,dataB.Size()+dataSB.Size());
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  fSB->SetParameters(result.GetParams()+1 );
  hSB->GetListOfFunctions()->Add(fSB);
  hSB->Draw(); 
  
  

}