rf804_mcstudy_constr.C: 'VALIDATION AND MC STUDIES' RooFit tutorial macro #804

/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #804
// 
// Using RooMCStudy on models with constrains
//
// 
// 07/2008 - Wouter Verkerke 
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;


void rf804_mcstudy_constr()
{
  // C r e a t e   m o d e l   w i t h   p a r a m e t e r   c o n s t r a i n t
  // ---------------------------------------------------------------------------

  // Observable
  RooRealVar x("x","x",-10,10) ;

  // Signal component
  RooRealVar m("m","m",0,-10,10) ;
  RooRealVar s("s","s",2,0.1,10) ;
  RooGaussian g("g","g",x,m,s) ;

  // Background component
  RooPolynomial p("p","p",x) ;

  // Composite model
  RooRealVar f("f","f",0.4,0.,1.) ;
  RooAddPdf sum("sum","sum",RooArgSet(g,p),f) ;

  // Construct constraint on parameter f
  RooGaussian fconstraint("fconstraint","fconstraint",f,RooConst(0.7),RooConst(0.1)) ;

  // Multiply constraint with p.d.f
  RooProdPdf sumc("sumc","sum with constraint",RooArgSet(sum,fconstraint)) ;



  // S e t u p   t o y   s t u d y   w i t h   m o d e l
  // ---------------------------------------------------

  // Perform toy study with internal constraint on f
  RooMCStudy mcs(sumc,x,Constrain(f),Silence(),Binned(),FitOptions(PrintLevel(-1))) ;

  // Run 500 toys of 2000 events.  
  // Before each toy is generated, a value for the f is sampled from the constraint pdf and 
  // that value is used for the generation of that toy.
  mcs.generateAndFit(500,2000) ;

  // Make plot of distribution of generated value of f parameter
  TH1* h_f_gen = mcs.fitParDataSet().createHistogram("f_gen",-40) ;

  // Make plot of distribution of fitted value of f parameter
  RooPlot* frame1  = mcs.plotParam(f,Bins(40)) ;
  frame1->SetTitle("Distribution of fitted f values") ;

  // Make plot of pull distribution on f
  RooPlot* frame2 = mcs.plotPull(f,Bins(40),FitGauss()) ;
  frame1->SetTitle("Distribution of f pull values") ;



  TCanvas* c = new TCanvas("rf804_mcstudy_constr","rf804_mcstudy_constr",1200,400) ;
  c->Divide(3) ;
  c->cd(1) ; h_f_gen->Draw() ;
  c->cd(2) ; frame1->Draw() ;
  c->cd(3) ; frame2->Draw() ;

}