Logo ROOT  
Reference Guide
rf801_mcstudy.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 ///
5 /// Validation and MC studies: toy Monte Carlo study that perform cycles of event generation and fitting
6 ///
7 /// \macro_image
8 /// \macro_output
9 /// \macro_code
10 ///
11 /// \date February 2018
12 /// \author Wouter Verkerke
13 
14 #include "RooRealVar.h"
15 #include "RooDataSet.h"
16 #include "RooGaussian.h"
17 #include "RooConstVar.h"
18 #include "RooChebychev.h"
19 #include "RooAddPdf.h"
20 #include "RooMCStudy.h"
21 #include "RooPlot.h"
22 #include "TCanvas.h"
23 #include "TAxis.h"
24 #include "TH2.h"
25 #include "RooFitResult.h"
26 #include "TStyle.h"
27 #include "TDirectory.h"
28 
29 using namespace RooFit;
30 
31 void rf801_mcstudy()
32 {
33  // C r e a t e m o d e l
34  // -----------------------
35 
36  // Declare observable x
37  RooRealVar x("x", "x", 0, 10);
38  x.setBins(40);
39 
40  // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
41  RooRealVar mean("mean", "mean of gaussians", 5, 0, 10);
42  RooRealVar sigma1("sigma1", "width of gaussians", 0.5);
43  RooRealVar sigma2("sigma2", "width of gaussians", 1);
44 
45  RooGaussian sig1("sig1", "Signal component 1", x, mean, sigma1);
46  RooGaussian sig2("sig2", "Signal component 2", x, mean, sigma2);
47 
48  // Build Chebychev polynomial p.d.f.
49  RooRealVar a0("a0", "a0", 0.5, 0., 1.);
50  RooRealVar a1("a1", "a1", -0.2, -1, 1.);
51  RooChebychev bkg("bkg", "Background", x, RooArgSet(a0, a1));
52 
53  // Sum the signal components into a composite signal p.d.f.
54  RooRealVar sig1frac("sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.);
55  RooAddPdf sig("sig", "Signal", RooArgList(sig1, sig2), sig1frac);
56 
57  // Sum the composite signal and background
58  RooRealVar nbkg("nbkg", "number of background events,", 150, 0, 1000);
59  RooRealVar nsig("nsig", "number of signal events", 150, 0, 1000);
60  RooAddPdf model("model", "g1+g2+a", RooArgList(bkg, sig), RooArgList(nbkg, nsig));
61 
62  // C r e a t e m a n a g e r
63  // ---------------------------
64 
65  // Instantiate RooMCStudy manager on model with x as observable and given choice of fit options
66  //
67  // The Silence() option kills all messages below the PROGRESS level, leaving only a single message
68  // per sample executed, and any error message that occur during fitting
69  //
70  // The Extended() option has two effects:
71  // 1) The extended ML term is included in the likelihood and
72  // 2) A poisson fluctuation is introduced on the number of generated events
73  //
74  // The FitOptions() given here are passed to the fitting stage of each toy experiment.
75  // If Save() is specified, the fit result of each experiment is saved by the manager
76  //
77  // A Binned() option is added in this example to bin the data between generation and fitting
78  // to speed up the study at the expense of some precision
79 
80  RooMCStudy *mcstudy =
82 
83  // G e n e r a t e a n d f i t e v e n t s
84  // ---------------------------------------------
85 
86  // Generate and fit 1000 samples of Poisson(nExpected) events
87  mcstudy->generateAndFit(1000);
88 
89  // E x p l o r e r e s u l t s o f s t u d y
90  // ------------------------------------------------
91 
92  // Make plots of the distributions of mean, the error on mean and the pull of mean
93  RooPlot *frame1 = mcstudy->plotParam(mean, Bins(40));
94  RooPlot *frame2 = mcstudy->plotError(mean, Bins(40));
95  RooPlot *frame3 = mcstudy->plotPull(mean, Bins(40), FitGauss(kTRUE));
96 
97  // Plot distribution of minimized likelihood
98  RooPlot *frame4 = mcstudy->plotNLL(Bins(40));
99 
100  // Make some histograms from the parameter dataset
101  TH1 *hh_cor_a0_s1f = mcstudy->fitParDataSet().createHistogram("hh", a1, YVar(sig1frac));
102  TH1 *hh_cor_a0_a1 = mcstudy->fitParDataSet().createHistogram("hh", a0, YVar(a1));
103 
104  // Access some of the saved fit results from individual toys
105  TH2 *corrHist000 = mcstudy->fitResult(0)->correlationHist("c000");
106  TH2 *corrHist127 = mcstudy->fitResult(127)->correlationHist("c127");
107  TH2 *corrHist953 = mcstudy->fitResult(953)->correlationHist("c953");
108 
109  // Draw all plots on a canvas
110  gStyle->SetOptStat(0);
111  TCanvas *c = new TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900);
112  c->Divide(3, 3);
113  c->cd(1);
114  gPad->SetLeftMargin(0.15);
115  frame1->GetYaxis()->SetTitleOffset(1.4);
116  frame1->Draw();
117  c->cd(2);
118  gPad->SetLeftMargin(0.15);
119  frame2->GetYaxis()->SetTitleOffset(1.4);
120  frame2->Draw();
121  c->cd(3);
122  gPad->SetLeftMargin(0.15);
123  frame3->GetYaxis()->SetTitleOffset(1.4);
124  frame3->Draw();
125  c->cd(4);
126  gPad->SetLeftMargin(0.15);
127  frame4->GetYaxis()->SetTitleOffset(1.4);
128  frame4->Draw();
129  c->cd(5);
130  gPad->SetLeftMargin(0.15);
131  hh_cor_a0_s1f->GetYaxis()->SetTitleOffset(1.4);
132  hh_cor_a0_s1f->Draw("box");
133  c->cd(6);
134  gPad->SetLeftMargin(0.15);
135  hh_cor_a0_a1->GetYaxis()->SetTitleOffset(1.4);
136  hh_cor_a0_a1->Draw("box");
137  c->cd(7);
138  gPad->SetLeftMargin(0.15);
139  corrHist000->GetYaxis()->SetTitleOffset(1.4);
140  corrHist000->Draw("colz");
141  c->cd(8);
142  gPad->SetLeftMargin(0.15);
143  corrHist127->GetYaxis()->SetTitleOffset(1.4);
144  corrHist127->Draw("colz");
145  c->cd(9);
146  gPad->SetLeftMargin(0.15);
147  corrHist953->GetYaxis()->SetTitleOffset(1.4);
148  corrHist953->Draw("colz");
149 
150  // Make RooMCStudy object available on command line after
151  // macro finishes
152  gDirectory->Add(mcstudy);
153 }
c
#define c(i)
Definition: RSha256.hxx:119
RooDataSet::createHistogram
TH2F * createHistogram(const RooAbsRealLValue &var1, const RooAbsRealLValue &var2, const char *cuts="", const char *name="hist") const
Create a TH2F histogram of the distribution of the specified variable using this dataset.
Definition: RooDataSet.cxx:1419
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooChebychev.h
RooAddPdf
Definition: RooAddPdf.h:32
TDirectory.h
RooFit::Bins
RooCmdArg Bins(Int_t nbin)
Definition: RooGlobalFunc.cxx:174
RooMCStudy::plotParam
RooPlot * plotParam(const RooRealVar &param, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Plot the distribution of the fitted value of the given parameter on a newly created frame.
Definition: RooMCStudy.cxx:1116
RooMCStudy.h
RooMCStudy::fitResult
const RooFitResult * fitResult(Int_t sampleNum) const
Return the RooFitResult of the fit with the given run number.
Definition: RooMCStudy.cxx:1019
RooFit::PrintEvalErrors
RooCmdArg PrintEvalErrors(Int_t numErrors)
Definition: RooGlobalFunc.cxx:205
RooChebychev
Definition: RooChebychev.h:25
RooMCStudy
Definition: RooMCStudy.h:32
RooMCStudy::plotPull
RooPlot * plotPull(const RooRealVar &param, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Plot the distribution of pull values for the specified parameter on a newly created frame.
Definition: RooMCStudy.cxx:1217
RooArgList
Definition: RooArgList.h:21
RooGaussian.h
TStyle.h
x
Double_t x[n]
Definition: legend1.C:17
RooGaussian
Definition: RooGaussian.h:25
RooAddPdf.h
TCanvas.h
RooFit::YVar
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
Definition: RooGlobalFunc.cxx:240
RooFitResult::correlationHist
TH2 * correlationHist(const char *name="correlation_matrix") const
Return TH2D of correlation matrix.
Definition: RooFitResult.cxx:1082
RooDataSet.h
RooMCStudy::plotNLL
RooPlot * plotNLL(const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Plot the distribution of the -log(L) values on a newly created frame.
Definition: RooMCStudy.cxx:1150
RooFit::FitOptions
RooCmdArg FitOptions(const char *opts)
Definition: RooGlobalFunc.cxx:184
RooFit::FitGauss
RooCmdArg FitGauss(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:282
gStyle
R__EXTERN TStyle * gStyle
Definition: TStyle.h:412
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:318
RooFit
Definition: RooCFunction1Binding.h:29
RooPlot.h
gDirectory
#define gDirectory
Definition: TDirectory.h:236
RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1256
TH2
Definition: TH2.h:30
RooPlot
Definition: RooPlot.h:44
RooMCStudy::fitParDataSet
const RooDataSet & fitParDataSet()
Return a RooDataSet containing the post-fit parameters of each toy cycle.
Definition: RooMCStudy.cxx:981
RooRealVar.h
TH2.h
RooFit::Binned
RooCmdArg Binned(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:268
TStyle::SetOptStat
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1592
RooFitResult.h
RooConstVar.h
TCanvas
Definition: TCanvas.h:23
RooMCStudy::plotError
RooPlot * plotError(const RooRealVar &param, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Plot the distribution of the fit errors for the specified parameter on a newly created frame.
Definition: RooMCStudy.cxx:1174
TAxis.h
rf801_mcstudy
Definition: rf801_mcstudy.py:1
TH1
Definition: TH1.h:57
gPad
#define gPad
Definition: TVirtualPad.h:287
RooMCStudy::generateAndFit
Bool_t generateAndFit(Int_t nSamples, Int_t nEvtPerSample=0, Bool_t keepGenData=kFALSE, const char *asciiFilePat=0)
Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
Definition: RooMCStudy.cxx:660
make_cnn_model.model
model
Definition: make_cnn_model.py:6
RooFit::Silence
RooCmdArg Silence(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:260
TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:293
RooRealVar
Definition: RooRealVar.h:35
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooArgSet
Definition: RooArgSet.h:28
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2997
RooFit::Extended
RooCmdArg Extended(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:155