//+ 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 <vector>

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

/*
// 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
};

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];
      for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ];

      double p2[5]; 
      for (int i = 0; i < 5; ++i) p2[i] = par[iparSB[i] ];

      return (*fChi2_1)(p1) + (*fChi2_2)(p2);
   } 

   const  ROOT::Math::IMultiGenFunction * fChi2_1;
   const  ROOT::Math::IMultiGenFunction * fChi2_2;
};
*/

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

  GlobalChi2(  const int /* length */,
	       ROOT::Math::IMultiGenFunction & farray) :
    fChi2_array (&farray)
  {}

  GlobalChi2(  TH1D *harray [2] ):
    hChi2_array (*harray)
  {}

  // 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];
    for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ];
    
    double p2[5];
    for (int i = 0; i < 5; ++i) p2[i] = par[iparSB[i] ];
    
    return (*fChi2_1)(p1) + (*fChi2_2)(p2);
  }

  void SetSpar(){
    for(int i = 1; i<6; i++)iparSB[i-1]=i;
  }

  void SetBpar(){
    iparB[0]=0;
    iparB[1]=2;
  }

  int *Spara(){
    return iparSB;
  }

  int *Bpara(){
    return iparB;
  }

  TH1D* getFirst(){
    return &hChi2_array[0];
  }

private:
  const  ROOT::Math::IMultiGenFunction * fChi2_1;
  const  ROOT::Math::IMultiGenFunction * fChi2_2;
  const  ROOT::Math::IMultiGenFunction * fChi2_array;
  TH1D * hChi2_array;
  int iparSB[5];
  int iparB[2];
};



int testcombinedFit_LL(void)
{ 



  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::DataOptions opt; 
  ROOT::Fit::DataRange rangeB; 
  // set the data range
  rangeB.SetRange(10,90);
  ROOT::Fit::BinData dataB(opt,rangeB); 
  ROOT::Fit::FillData(dataB, hB);

  ROOT::Fit::DataRange rangeSB; 
  rangeSB.SetRange(10,50);
  ROOT::Fit::BinData dataSB(opt,rangeSB); 
  ROOT::Fit::FillData(dataSB, hSB);

  ROOT::Fit::PoissonLLFunction ll_B(dataB, wfB);
  ROOT::Fit::PoissonLLFunction ll_SB(dataSB, wfSB);

  ROOT::Fit::PoissonLLFunction *test[2];
  test[0] = new ROOT::Fit::PoissonLLFunction(dataB, wfB);
  test[1] = new ROOT::Fit::PoissonLLFunction(dataSB, wfSB);

    //  std::vector<ROOT::Fit::PoissonLLFunction> test;
  //  test.resize(2);
  //  test.push_back(ll_B);
  //  test.push_back(ll_SB);

  GlobalChi2 globalLL(ll_B, ll_SB);
  globalLL.SetSpar();
  globalLL.SetBpar();

  GlobalChi2 testLL(2, **test);

  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 if it is a chi2 fit)
  fitter.FitFCN(6,globalLL,0,dataB.Size()+dataSB.Size(),false);
  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,2);
  c1->cd(1);
  gStyle->SetOptFit(1111);

  fB->SetFitResult( result, globalLL.Bpara());
  fB->SetRange(rangeB().first, rangeB().second);   
  fB->SetLineColor(kBlue);
  hB->GetListOfFunctions()->Add(fB);
  hB->Draw(); 

  c1->cd(2);
  fSB->SetFitResult( result, globalLL.Spara());
  fSB->SetRange(rangeSB().first, rangeSB().second);   
  fSB->SetLineColor(kRed);
  hSB->GetListOfFunctions()->Add(fSB);
  hSB->Draw(); 

  TH1D *htest[2];
  htest[0] = new TH1D("blub0", "", 1000, 0, 100);
  htest[0]->FillRandom("fS",100);
  htest[1] = new TH1D("blub1", "", 100, 0, 100);
  htest[1]->FillRandom("fS",100);

  
  GlobalChi2 histtest (htest);

  TH1D* blub = histtest.getFirst();

  TCanvas *c3 = new TCanvas("c3","c3",1);
  c3 = c3; // get rid of compiler's warning: unused variable ‘c3’
  blub->Draw();

#if 1 /* 0 or 1 */
  TFile *rfile = new TFile("blub.root", "recreate");
#else /* 0 or 1 */
  TFile *rfile = TFile::Open("blub.root", "recreate");
#endif /* 0 or 1 */
  blub->Write();
  rfile->Close();
  delete rfile;

  return 0;
}

#if !defined(__CINT__) && !defined(__ACLIC__)
int main(int argc, char** argv)
{
  // TApplication app("app", &argc, argv);
  
  return testcombinedFit_LL();
}
#endif