#ifdef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooAbsPdf.h"
#include "RooMinimizer.h"
#include "RooAbsReal.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooExtendPdf.h"
#include "RooAddPdf.h"
#include "TRandom3.h"

#include <vector>

std::vector<RooDataSet*> makePseudoSamples()
{
  RooRealVar x("x","x",-10,10);
  RooRealVar m("m","m",0,-10,10);
  RooRealVar s("s","s",1,0  ,10);
  RooGaussian gaus("gaus","gaus",x,m,s);
  RooRealVar aCheb("aCheb","aCheb",-0.05,-10,10);
  RooChebychev cheb("cheb","cheb",x,aCheb);
  
  RooRealVar nGaus("nGaus","nGaus",25,-1000,1000);
  RooRealVar nCheb("nCheb","nCheb",30,-1000,1000);
  
  RooAddPdf combPDF("combPDF","combPDF",RooArgList(gaus,cheb),RooArgList(nGaus,nCheb));

  TRandom3 *rand = new TRandom3();
  std::vector<RooDataSet*> retSets;
  for(int iSet = 0; iSet < 10; iSet++)
    {
      char setName[100];
      sprintf(setName,"toySet%d",iSet);
      RooDataSet *tempSet = combPDF.generate(x,rand->Poisson(55));
      retSets.push_back((RooDataSet*)tempSet->Clone(setName));
      delete tempSet;
    }
  return retSets;
}

int fitSample(RooAbsPdf *myPDF, RooDataSet *mySet)
{
  RooAbsReal *nll = myPDF->createNLL(*mySet,RooFit::Extended(true),RooFit::Range("fullRange"));
  RooMinimizer min(*nll);
  min.setPrintLevel(3);
  int migradStat = min.migrad();
  int hesseStat  = min.hesse();
  delete nll;
  if(migradStat != 0 || hesseStat != 0) printf("MIGRAD: %d HESSE: %d\n",migradStat,hesseStat);
  return migradStat + hesseStat;
}

int testIterativeFit()
{
  RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
  RooMsgService::instance().setSilentMode(true);
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  RooMsgService::instance().setSilentMode(true);
  std::vector<RooDataSet*> toySamples = makePseudoSamples();
  
  RooRealVar *x        = new RooRealVar("x"     ,"x"    ,-10,10);
  RooRealVar *mean     = new RooRealVar("mean"  ,"mean" ,0  ,-10  ,10);
  RooRealVar *sigma    = new RooRealVar("sigma" ,"sigma",1  ,0    ,10);
  RooRealVar *nGaus    = new RooRealVar("nGaus" ,"nGaus",25,-1000,1000);
  RooGaussian *myGaus  = new RooGaussian("gaus" ,"gaus" ,*x,*mean,*sigma);
  RooRealVar *aCheb    = new RooRealVar("aCheb" ,"aCheb",-0.05,-10,10);
  RooRealVar *nCheb    = new RooRealVar("nCheb" ,"nCheb",30,-1000,1000);
  RooChebychev *myCheb = new RooChebychev("cheb","cheb" ,*x,*aCheb);

  RooAddPdf *combPDF   = new RooAddPdf("combPDF","combPDF",RooArgList(*myGaus,*myCheb),RooArgSet(*nGaus,*nCheb));

  int oneFit = 9;
  oneFit = 0; //if you comment this line all 10 fits are going to be run otherwise only the last one to show the standalone part
  for(int iToy = oneFit; iToy < (int)toySamples.size(); iToy++)
    {
      RooDataSet *curSet   = toySamples[iToy];
      const RooArgSet *obs = curSet->get();
      x                    = (RooRealVar*)obs->find("x");
      x->setRange("fullRange",-10,10);
      printf("TOY %d FITTED!\n",iToy);
      int fitStatus        = fitSample((RooAbsPdf*)combPDF,curSet);
      //printf("Set: %d status = %d\n",iToy,fitStatus);
      mean ->setVal(0)    ; mean ->setError(0);
      sigma->setVal(1)    ; sigma->setError(0);
      aCheb->setVal(-0.05); aCheb->setError(0);
      nGaus->setVal(100)  ; nGaus->setError(0);
      nCheb->setVal(200)  ; nCheb->setError(0);
      delete toySamples[iToy];
    }

  return -1;
}