/////////////////////////////////////////////////////////////////////////
//
// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #504
// 
// Using RooSimWSTool to construct a simultaneous p.d.f that is built
// of variations of an input p.d.f
//
//
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooCategory.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooPolynomial.h"
#include "RooSimultaneous.h"
#include "RooAddPdf.h"
#include "RooWorkspace.h"
#include "RooSimWSTool.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TFile.h"
#include "TH1.h"
using namespace RooFit ;


void rf504_simwstool()
{
  // C r e a t e   m a s t e r   p d f 
  // ---------------------------------
  
  // Construct gauss(x,m,s)
  RooRealVar x("x","x",-10,10) ;
  RooRealVar m("m","m",0,-10,10) ;
  RooRealVar s("s","s",1,-10,10) ;
  RooGaussian gauss("g","g",x,m,s) ;

  // Construct poly(x,p0)
  RooRealVar p0("p0","p0",0.01,0.,1.) ;
  RooPolynomial poly("p","p",x,p0) ;		 

  // Construct model = f*gauss(x) + (1-f)*poly(x)
  RooRealVar f("f","f",0.5,0.,1.) ;
  RooAddPdf model("model","model",RooArgSet(gauss,poly),f) ;



  // C r e a t e   c a t e g o r y   o b s e r v a b l e s   f o r   s p l i t t i n g
  // ----------------------------------------------------------------------------------

  // Define two categories that can be used for splitting
  RooCategory c("c","c") ;
  c.defineType("run1") ;
  c.defineType("run2") ;

  RooCategory d("d","d") ;
  d.defineType("foo") ;
  d.defineType("bar") ;



  // S e t u p   S i m W S T o o l 
  // -----------------------------

  // Import ingredients in a workspace
  RooWorkspace w("w","w") ;
  w.import(RooArgSet(model,c,d)) ;

  // Make Sim builder tool
  RooSimWSTool sct(w) ;


  // B u i l d   a   s i m u l t a n e o u s   m o d e l   w i t h   o n e   s p l i t
  // ---------------------------------------------------------------------------------

  // Construct a simultaneous p.d.f with the following form
  //
  // model_run1(x) = f*gauss_run1(x,m_run1,s) + (1-f)*poly
  // model_run2(x) = f*gauss_run2(x,m_run2,s) + (1-f)*poly
  // simpdf(x,c) = model_run1(x) if c=="run1"
  //             = model_run2(x) if c=="run2"
  //
  // Returned p.d.f is owned by the workspace
  RooSimultaneous* model_sim = sct.build("model_sim","model",SplitParam("m","c")) ;

  // Print tree structure of model
  model_sim->Print("t") ;

  // Adjust model_sim parameters in workspace
  w.var("m_run1")->setVal(-3) ;
  w.var("m_run2")->setVal(+3) ;

  // Print contents of workspace
  w.Print("v") ;



  // B u i l d   a   s i m u l t a n e o u s   m o d e l   w i t h   p r o d u c t   s p l i t
  // -----------------------------------------------------------------------------------------

  // Build another simultaneous p.d.f using a composite split in states c X d
  RooSimultaneous* model_sim2 = sct.build("model_sim2","model",SplitParam("p0","c,d")) ;

  // Print tree structure of this model
  model_sim2->Print("t") ;

}