ROOT logo

From $ROOTSYS/tutorials/roofit/rf801_mcstudy.C

/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #801
// 
// A Toy Monte Carlo study that perform cycles of
// event generation and fittting
//
// 
/////////////////////////////////////////////////////////////////////////

#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 "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH2.h"
#include "RooFitResult.h"
#include "TStyle.h"
#include "TDirectory.h"

using namespace RooFit ;


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

  // Declare observable x
  RooRealVar x("x","x",0,10) ;
  x.setBins(40) ;

  // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their paramaters
  RooRealVar mean("mean","mean of gaussians",5,0,10) ;
  RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
  RooRealVar sigma2("sigma2","width of gaussians",1) ;

  RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;  
  RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;  
  
  // Build Chebychev polynomial p.d.f.  
  RooRealVar a0("a0","a0",0.5,0.,1.) ;
  RooRealVar a1("a1","a1",-0.2,-1,1.) ;
  RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;

  // Sum the signal components into a composite signal p.d.f.
  RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
  RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;

  // Sum the composite signal and background 
  RooRealVar nbkg("nbkg","number of background events,",150,0,1000) ;
  RooRealVar nsig("nsig","number of signal events",150,0,1000) ;
  RooAddPdf  model("model","g1+g2+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ;



  // C r e a t e   m a n a g e r
  // ---------------------------

  // Instantiate RooMCStudy manager on model with x as observable and given choice of fit options
  //
  // The Silence() option kills all messages below the PROGRESS level, leaving only a single message
  // per sample executed, and any error message that occur during fitting
  //
  // The Extended() option has two effects: 
  //    1) The extended ML term is included in the likelihood and 
  //    2) A poisson fluctuation is introduced on the number of generated events 
  //
  // The FitOptions() given here are passed to the fitting stage of each toy experiment.
  // If Save() is specified, the fit result of each experiment is saved by the manager  
  //
  // A Binned() option is added in this example to bin the data between generation and fitting
  // to speed up the study at the expemse of some precision

  RooMCStudy* mcstudy = new RooMCStudy(model,x,Binned(kTRUE),Silence(),Extended(),
				       FitOptions(Save(kTRUE),PrintEvalErrors(0))) ;
  

  // G e n e r a t e   a n d   f i t   e v e n t s
  // ---------------------------------------------

  // Generate and fit 1000 samples of Poisson(nExpected) events
  mcstudy->generateAndFit(1000) ;



  // E x p l o r e   r e s u l t s   o f   s t u d y 
  // ------------------------------------------------

  // Make plots of the distributions of mean, the error on mean and the pull of mean
  RooPlot* frame1 = mcstudy->plotParam(mean,Bins(40)) ;
  RooPlot* frame2 = mcstudy->plotError(mean,Bins(40)) ;
  RooPlot* frame3 = mcstudy->plotPull(mean,Bins(40),FitGauss(kTRUE)) ;

  // Plot distribution of minimized likelihood
  RooPlot* frame4 = mcstudy->plotNLL(Bins(40)) ;

  // Make some histograms from the parameter dataset
  TH1* hh_cor_a0_s1f = mcstudy->fitParDataSet().createHistogram("hh",a1,YVar(sig1frac)) ;
  TH1* hh_cor_a0_a1  = mcstudy->fitParDataSet().createHistogram("hh",a0,YVar(a1)) ;

  // Access some of the saved fit results from individual toys
  TH2* corrHist000 = mcstudy->fitResult(0)->correlationHist("c000") ;
  TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ;
  TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ;



  // Draw all plots on a canvas
  gStyle->SetPalette(1) ;
  gStyle->SetOptStat(0) ;
  TCanvas* c = new TCanvas("rf801_mcstudy","rf801_mcstudy",900,900) ;
  c->Divide(3,3) ;
  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() ;
  c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3->Draw() ;
  c->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.4) ; frame4->Draw() ;
  c->cd(5) ; gPad->SetLeftMargin(0.15) ; hh_cor_a0_s1f->GetYaxis()->SetTitleOffset(1.4) ; hh_cor_a0_s1f->Draw("box") ;
  c->cd(6) ; gPad->SetLeftMargin(0.15) ; hh_cor_a0_a1->GetYaxis()->SetTitleOffset(1.4) ; hh_cor_a0_a1->Draw("box") ;
  c->cd(7) ; gPad->SetLeftMargin(0.15) ; corrHist000->GetYaxis()->SetTitleOffset(1.4) ; corrHist000->Draw("colz") ;
  c->cd(8) ; gPad->SetLeftMargin(0.15) ; corrHist127->GetYaxis()->SetTitleOffset(1.4) ; corrHist127->Draw("colz") ;
  c->cd(9) ; gPad->SetLeftMargin(0.15) ; corrHist953->GetYaxis()->SetTitleOffset(1.4) ; corrHist953->Draw("colz") ;

  // Make RooMCStudy object available on command line after
  // macro finishes
  gDirectory->Add(mcstudy) ;
}
 rf801_mcstudy.C:1
 rf801_mcstudy.C:2
 rf801_mcstudy.C:3
 rf801_mcstudy.C:4
 rf801_mcstudy.C:5
 rf801_mcstudy.C:6
 rf801_mcstudy.C:7
 rf801_mcstudy.C:8
 rf801_mcstudy.C:9
 rf801_mcstudy.C:10
 rf801_mcstudy.C:11
 rf801_mcstudy.C:12
 rf801_mcstudy.C:13
 rf801_mcstudy.C:14
 rf801_mcstudy.C:15
 rf801_mcstudy.C:16
 rf801_mcstudy.C:17
 rf801_mcstudy.C:18
 rf801_mcstudy.C:19
 rf801_mcstudy.C:20
 rf801_mcstudy.C:21
 rf801_mcstudy.C:22
 rf801_mcstudy.C:23
 rf801_mcstudy.C:24
 rf801_mcstudy.C:25
 rf801_mcstudy.C:26
 rf801_mcstudy.C:27
 rf801_mcstudy.C:28
 rf801_mcstudy.C:29
 rf801_mcstudy.C:30
 rf801_mcstudy.C:31
 rf801_mcstudy.C:32
 rf801_mcstudy.C:33
 rf801_mcstudy.C:34
 rf801_mcstudy.C:35
 rf801_mcstudy.C:36
 rf801_mcstudy.C:37
 rf801_mcstudy.C:38
 rf801_mcstudy.C:39
 rf801_mcstudy.C:40
 rf801_mcstudy.C:41
 rf801_mcstudy.C:42
 rf801_mcstudy.C:43
 rf801_mcstudy.C:44
 rf801_mcstudy.C:45
 rf801_mcstudy.C:46
 rf801_mcstudy.C:47
 rf801_mcstudy.C:48
 rf801_mcstudy.C:49
 rf801_mcstudy.C:50
 rf801_mcstudy.C:51
 rf801_mcstudy.C:52
 rf801_mcstudy.C:53
 rf801_mcstudy.C:54
 rf801_mcstudy.C:55
 rf801_mcstudy.C:56
 rf801_mcstudy.C:57
 rf801_mcstudy.C:58
 rf801_mcstudy.C:59
 rf801_mcstudy.C:60
 rf801_mcstudy.C:61
 rf801_mcstudy.C:62
 rf801_mcstudy.C:63
 rf801_mcstudy.C:64
 rf801_mcstudy.C:65
 rf801_mcstudy.C:66
 rf801_mcstudy.C:67
 rf801_mcstudy.C:68
 rf801_mcstudy.C:69
 rf801_mcstudy.C:70
 rf801_mcstudy.C:71
 rf801_mcstudy.C:72
 rf801_mcstudy.C:73
 rf801_mcstudy.C:74
 rf801_mcstudy.C:75
 rf801_mcstudy.C:76
 rf801_mcstudy.C:77
 rf801_mcstudy.C:78
 rf801_mcstudy.C:79
 rf801_mcstudy.C:80
 rf801_mcstudy.C:81
 rf801_mcstudy.C:82
 rf801_mcstudy.C:83
 rf801_mcstudy.C:84
 rf801_mcstudy.C:85
 rf801_mcstudy.C:86
 rf801_mcstudy.C:87
 rf801_mcstudy.C:88
 rf801_mcstudy.C:89
 rf801_mcstudy.C:90
 rf801_mcstudy.C:91
 rf801_mcstudy.C:92
 rf801_mcstudy.C:93
 rf801_mcstudy.C:94
 rf801_mcstudy.C:95
 rf801_mcstudy.C:96
 rf801_mcstudy.C:97
 rf801_mcstudy.C:98
 rf801_mcstudy.C:99
 rf801_mcstudy.C:100
 rf801_mcstudy.C:101
 rf801_mcstudy.C:102
 rf801_mcstudy.C:103
 rf801_mcstudy.C:104
 rf801_mcstudy.C:105
 rf801_mcstudy.C:106
 rf801_mcstudy.C:107
 rf801_mcstudy.C:108
 rf801_mcstudy.C:109
 rf801_mcstudy.C:110
 rf801_mcstudy.C:111
 rf801_mcstudy.C:112
 rf801_mcstudy.C:113
 rf801_mcstudy.C:114
 rf801_mcstudy.C:115
 rf801_mcstudy.C:116
 rf801_mcstudy.C:117
 rf801_mcstudy.C:118
 rf801_mcstudy.C:119
 rf801_mcstudy.C:120
 rf801_mcstudy.C:121
 rf801_mcstudy.C:122
 rf801_mcstudy.C:123
 rf801_mcstudy.C:124
 rf801_mcstudy.C:125
 rf801_mcstudy.C:126
 rf801_mcstudy.C:127
 rf801_mcstudy.C:128
 rf801_mcstudy.C:129
 rf801_mcstudy.C:130
 rf801_mcstudy.C:131
 rf801_mcstudy.C:132
 rf801_mcstudy.C:133
 rf801_mcstudy.C:134
 rf801_mcstudy.C:135
 rf801_mcstudy.C:136