//+ 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 "TH1.h"
#include "TF1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.h"
#include <cmath> 
#include "TGraph.h"
#include "TRandom.h"

// definition of shared parameter
// background function 
int iparB1[2] = { 0,      // exp amplitude in B histo
                 2    // exp common parameter 
};

// signal + background function 
int iparB2[2] = { 1, // exp amplitude in S+B histo
                  2, // exp common parameter
};

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[iparB1[i] ];

      double p2[2]; 
      for (int i = 0; i < 2; ++i) p2[i] = par[iparB2[i] ];

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

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

double myfunc1(double *x, double *par){
      float xx =x[0];
      double f = exp(par[0]*xx)*sin(par[1]*xx);
      return f;
   }
   
double myfunc2(double *x, double *par){
      float xx =x[0];
      double f = exp(par[0]*xx)*sin(par[1]*xx);
      return f;
   }
   
double myfunction1(double *x, double *par){ 
      double constant = 1.0;
      float xx =x[0];
      double f = exp(par[0]*xx)*sin(par[1]*xx);
      return f;
   }
   
double myfunction2(double *x, double *par){
      float constant = 1.0;
      float xx =x[0];
      double f = exp(par[0]*xx)*sin(par[1]*xx);
      return f;
   }

void combinedFit(){ 
  int nd = 100;
  double rnd1,x, rnd2,xp[nd],y1[nd],y2[nd];
  double e = 0.2;
  TRandom r1,r2;
  x = (double) 100.0/nd;
  
  //TH1D * hB1 = new TH1D("hB1","histo B1",100,0,100);
  //TH1D * hB2 = new TH1D("hB2","histo B2",100, 0,100);

  TF1 * fB1 = new TF1("fB1",myfunc1,0,100,2);
  TF1 * fB2 = new TF1("fB2",myfunc2,0,100,2);
  
  fB1->SetParameters(-0.05,1);
  fB2->SetParameters(-0.15,1);
  
  //hB1->FillRandom("fB1",2000);
  //hB2->FillRandom("fB1",2000);
  for (Int_t i=0; i<nd; i++) {
      rnd1 = r1.Uniform(-e,e); // Generate a random number in [-e,e]
      rnd2 = r2.Uniform(-e,e); // Generate a random number in [-e,e]
      y1[i] = fB1->Eval(x*i)*(1+rnd1);
      y2[i] = fB2->Eval(x*i)*(1+rnd2);
      xp[i] = x*i;
      //cout<< rnd2<<" "<<xp[i] <<" " << y2[i]<<endl;
  }
  
  TGraph *gb1 = new TGraph(nd,xp,y1);
  TGraph *gb2 = new TGraph(nd,xp,y2);
  
  TF1 * ffitB1 = new TF1("ffitB1",myfunction1,0,100,2);
  TF1 * ffitB2 = new TF1("ffitB2",myfunction2,0,100,2);
  
  ROOT::Math::WrappedMultiTF1 wfB1(*ffitB1,1);
  ROOT::Math::WrappedMultiTF1 wfB2(*ffitB2,1);

  ROOT::Fit::DataOptions opt; 
  ROOT::Fit::DataRange rangeB1; 
  ROOT::Fit::DataRange rangeB2; 
 
  // set the data range
  rangeB1.SetRange(0,100);
  rangeB2.SetRange(0,100);
 
  ROOT::Fit::BinData dataB1(opt,rangeB1); 
  ROOT::Fit::BinData dataB2(opt,rangeB2); 
 
  ROOT::Fit::FillData(dataB1, gb1);
  ROOT::Fit::FillData(dataB2, gb2);

 
  ROOT::Fit::Chi2Function chi2_B1(dataB1, wfB1);
  ROOT::Fit::Chi2Function chi2_B2(dataB2, wfB2);

  GlobalChi2 globalChi2(chi2_B1, chi2_B2);

  ROOT::Fit::Fitter fitter;

  const int Npar = 3; 
  double par0[Npar] = {0.1,0.1,5.0};

  // create before the parameter settings in order to fix or set range on them
  fitter.Config().SetParamsSettings(3,par0);
  fitter.Config().MinimizerOptions().SetPrintLevel(1);
  fitter.Config().SetMinimizer("Minuit","Minimize"); 

  // fit FCN function directly 
  // (specify optionally data size and flag to indicate that is a chi2 fit)
  fitter.FitFCN(3,globalChi2,0,dataB1.Size()+dataB2.Size(),true);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two graphs",
                             10,10,700,700);
  c1->Divide(1,2);
  c1->cd(1);
  gStyle->SetOptFit(1111);

  ffitB1->SetFitResult(result, iparB1);
  ffitB1->SetRange(rangeB1().first, rangeB1().second);   
  ffitB1->SetLineColor(kBlue);
  gb1->GetListOfFunctions()->Add(ffitB1);
  gb1->Draw("alp"); 
// 
  c1->cd(2);
  ffitB2->SetFitResult( result, iparB2);
  ffitB2->SetRange(rangeB2().first, rangeB2().second);   
  ffitB2->SetLineColor(kRed);
  gb2->GetListOfFunctions()->Add(ffitB2);
  gb2->Draw("alp"); 

//TCanvas * c1 = new TCanvas("c1");
//   c1->Divide(1,2);
//   c1->cd(1);
//   
//   gb1->Draw("alp");
//   c1->cd(2);
//   
//   gb2->Draw("alp");

}