//+ 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 "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.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(  vector< ROOT::Math::IMultiGenFunction * > & fVec) : 
      fChi2_Vec(&fVec) {}

   // parameter vector is first background (in common 1 and 2) 
   // and then is signal (only in 2)
   double operator() (const double *par) const {
      vector< vector<double> > pVec;
      vector< double > dummyVec;
      for (int i = 0; i < 2; ++i) dummyVec.push_back(par[iparB[i] ]);
      pVec.push_back(dummyVec);

      dummyVec.clear();
      for (int i = 0; i < 5; ++i) dummyVec.push_back(par[iparSB[i] ]);
      pVec.push_back(dummyVec);
      
      double val = 0;
      for( size_t i = 0; i < (*fChi2_Vec).size(); i++ ){
        val += (*fChi2_Vec->at(i))(&(pVec[i][0]));
      }

      return val;
   } 

   vector< ROOT::Math::IMultiGenFunction * > * fChi2_Vec;
};

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::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);
  
  vector< ROOT::Fit::Chi2Function * > chi2_vec;
  chi2_vec.push_back( new ROOT::Fit::Chi2Function(dataB, wfB) );
  chi2_vec.push_back( new ROOT::Fit::Chi2Function(dataSB, wfSB) );

  GlobalChi2 globalChi2(chi2_vec);
}