/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #803
// 
// RooMCStudy: Using the randomizer and profile likelihood add-on models
//
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooRandomizeParamMCSModule.h"
#include "RooDLLSignificanceMCSModule.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TDirectory.h"

using namespace RooFit ;


void MCStudyEg()
{
  // C r e a t e   m o d e l 
  // -----------------------

  // Simulation of signal and background of top quark decaying into
  // 3 jets with background

  // Observable
  RooRealVar mjjj("mjjj","m(3jet) (GeV)",100,85.,350.) ;

  // Signal component (Gaussian)
  RooRealVar mtop("mtop","m(top)",162) ;
  RooRealVar wtop("wtop","m(top) resolution",15.2) ;
  RooGaussian sig("sig","top signal",mjjj,mtop,wtop) ;

  // Background component (Chebychev)
  RooRealVar c0("c0","Chebychev coefficient 0",-0.846,-1.,1.) ;
  RooRealVar c1("c1","Chebychev coefficient 1", 0.112,-1.,1.) ;
  RooRealVar c2("c2","Chebychev coefficient 2", 0.076,-1.,1.) ;
  RooChebychev bkg("bkg","combinatorial background",mjjj,RooArgList(c0,c1,c2)) ;

  // Composite model
  RooRealVar nsig("nsig","number of signal events",53,0,1e3) ;
  RooRealVar nbkg("nbkg","number of background events",103,0,5e3) ;
  RooAddPdf model("model","model",RooArgList(sig,bkg),RooArgList(nsig,nbkg)) ;

  RooMCStudy* mcs = new RooMCStudy(model,mjjj,Binned(),Silence(),Extended(kTRUE),
				   FitOptions(Extended(kTRUE),PrintEvalErrors(-1))) ;

  
  RooRandomizeParamMCSModule randModule ;
  //randModule.sampleSumUniform(nsig,50,500) ;
  randModule.sampleGaussian(nsig,120,1);
  mcs->addModule(randModule) ;   

  mcs->generateAndFit(200) ;

  TCanvas *c = new TCanvas("NFitted_sig", "canvas",640,480);
  RooPlot* frame1 = nsig.frame(Name("NEvents"),Title("Number of signal events"));
  mcs.plotParamOn(frame1);
  c->cd(1); frame1->Draw();
  c->SaveAs("NFitted_sig.eps");
  c->Write();
 
  gDirectory->Add(mcs) ;

}