#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit ;


void SimultaneousTest(Bool_t manuFit = false)
{
  // C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e
  // -------------------------------------------------------------
  
  // Create observables
  RooRealVar x("x","x",-8,8) ;

  // Construct signal pdf
  RooRealVar mean("mean","mean",0,-8,8) ;
  RooRealVar sigma("sigma","sigma",0.3,0.1,10) ;
  RooGaussian gx("gx","gx",x,mean,sigma) ;

  // Construct background pdf
  RooRealVar a0("a0","a0",-0.1,-1,1) ;
  RooRealVar a1("a1","a1",0.004,-1,1) ;
  RooChebychev px("px","px",x,RooArgSet(a0,a1)) ;

  // Construct composite pdf
  RooRealVar S("S","S",9.,0.,100000000.) ;
  RooRealVar B("B","B",1.,0.,100000000.) ;
  RooAddPdf model("model","model",RooArgList(gx,px),RooArgList(S,B)) ;



  // C r e a t e   m o d e l   f o r   c o n t r o l   s a m p l e
  // --------------------------------------------------------------

  // Construct the composite model
  RooRealVar S_ctl("S_ctl","S_ctl",9.,0.,100000000.) ;
  RooRealVar B_ctl("B_ctl","B_ctl",1.,0.,100000000.) ;
  RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx,px),RooArgList(S_ctl,B_ctl)) ;



  // G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s 
  // ---------------------------------------------------------------

  // Generate events from model
  RooDataSet *data = model.generate(RooArgSet(x),1000) ;
  RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x),10000) ;



  // C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s 
  // ---------------------------------------------------------------------------

  // Define category to distinguish physics and control samples events
  RooCategory sample("sample","sample") ;
  sample.defineType("physics") ;
  sample.defineType("control") ;

  // Construct combined dataset in (x,sample)
  RooDataSet combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ;



  // C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )
  // -----------------------------------------------------------------------------------

  // Construct a simultaneous pdf using category sample as index
  RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;

  // Associate model with the physics state and model_ctl with the control state
  simPdf.addPdf(model,"physics") ;
  simPdf.addPdf(model_ctl,"control") ;
  
  x.setRange("fitrange_physics", -8., 8.);
  x.setRange("fitrange_control", -8., 3.); 


  
  // P e r f o r m   a   s i m u l t a n e o u s   f i t
  // ---------------------------------------------------
  
  // Perform simultaneous fit of model to data and model_ctl to data_ctl
  RooFitResult *fitResult = simPdf.fitTo(combData, Range("fitrange"), SplitRange(kTRUE), Extended(kTRUE), Verbose(kTRUE), Save(true)) ;
  cout << endl
       << "===================================" << endl
       << "=" << endl
       << "= Status  = " << fitResult->status()  << endl
       << "= Edm     = " << fitResult->edm()     << endl
       << "= CovQual = " << fitResult->covQual() << endl
       << "=" << endl
       << "===================================" << endl
       << endl;
  


  // P l o t   m o d e l   s l i c e s   o n   d a t a    s l i c e s 
  // ----------------------------------------------------------------

  // Make a frame for the physics sample
  RooPlot* frame1 = x.frame(Title("Physics sample")) ;

  // Plot all data tagged as physics sample
  combData.plotOn(frame1,Cut("sample==sample::physics")) ;

  // Plot "physics" slice of simultaneous pdf. 
  simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ;
  simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;

  // The same plot for the control sample slice
  // RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ;
  RooPlot* frame2 = x.frame(Title("Control sample")) ;
  combData.plotOn(frame2,Cut("sample==sample::control")) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;



  TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf",800,400) ;
  c->Divide(2) ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ;


}