ROOT   Reference Guide
Searching...
No Matches
rf801_mcstudy.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4/// Validation and MC studies: toy Monte Carlo study that perform cycles of event generation and fitting
5///
6/// \macro_image
7/// \macro_code
8/// \macro_output
9///
10/// \date February 2018
11/// \author Wouter Verkerke
12
13#include "RooRealVar.h"
14#include "RooDataSet.h"
15#include "RooGaussian.h"
16#include "RooChebychev.h"
18#include "RooMCStudy.h"
19#include "RooPlot.h"
20#include "TCanvas.h"
21#include "TAxis.h"
22#include "TH2.h"
23#include "RooFitResult.h"
24#include "TStyle.h"
25#include "TDirectory.h"
26
27using namespace RooFit;
28
29void rf801_mcstudy()
30{
31 // C r e a t e m o d e l
32 // -----------------------
33
34 // Declare observable x
35 RooRealVar x("x", "x", 0, 10);
36 x.setBins(40);
37
38 // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
39 RooRealVar mean("mean", "mean of gaussians", 5, 0, 10);
40 RooRealVar sigma1("sigma1", "width of gaussians", 0.5);
41 RooRealVar sigma2("sigma2", "width of gaussians", 1);
42
43 RooGaussian sig1("sig1", "Signal component 1", x, mean, sigma1);
44 RooGaussian sig2("sig2", "Signal component 2", x, mean, sigma2);
45
46 // Build Chebychev polynomial pdf
47 RooRealVar a0("a0", "a0", 0.5, 0., 1.);
48 RooRealVar a1("a1", "a1", -0.2, -1, 1.);
49 RooChebychev bkg("bkg", "Background", x, RooArgSet(a0, a1));
50
51 // Sum the signal components into a composite signal pdf
52 RooRealVar sig1frac("sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.);
53 RooAddPdf sig("sig", "Signal", RooArgList(sig1, sig2), sig1frac);
54
55 // Sum the composite signal and background
56 RooRealVar nbkg("nbkg", "number of background events,", 150, 0, 1000);
57 RooRealVar nsig("nsig", "number of signal events", 150, 0, 1000);
58 RooAddPdf model("model", "g1+g2+a", RooArgList(bkg, sig), RooArgList(nbkg, nsig));
59
60 // C r e a t e m a n a g e r
61 // ---------------------------
62
63 // Instantiate RooMCStudy manager on model with x as observable and given choice of fit options
64 //
65 // The Silence() option kills all messages below the PROGRESS level, leaving only a single message
66 // per sample executed, and any error message that occur during fitting
67 //
68 // The Extended() option has two effects:
69 // 1) The extended ML term is included in the likelihood and
70 // 2) A poisson fluctuation is introduced on the number of generated events
71 //
72 // The FitOptions() given here are passed to the fitting stage of each toy experiment.
73 // If Save() is specified, the fit result of each experiment is saved by the manager
74 //
75 // A Binned() option is added in this example to bin the data between generation and fitting
76 // to speed up the study at the expense of some precision
77
78 RooMCStudy *mcstudy =
79 new RooMCStudy(model, x, Binned(true), Silence(), Extended(), FitOptions(Save(true), PrintEvalErrors(0)));
80
81 // G e n e r a t e a n d f i t e v e n t s
82 // ---------------------------------------------
83
84 // Generate and fit 1000 samples of Poisson(nExpected) events
85 mcstudy->generateAndFit(1000);
86
87 // E x p l o r e r e s u l t s o f s t u d y
88 // ------------------------------------------------
89
90 // Make plots of the distributions of mean, the error on mean and the pull of mean
91 RooPlot *frame1 = mcstudy->plotParam(mean, Bins(40));
92 RooPlot *frame2 = mcstudy->plotError(mean, Bins(40));
93 RooPlot *frame3 = mcstudy->plotPull(mean, Bins(40), FitGauss(true));
94
95 // Plot distribution of minimized likelihood
96 RooPlot *frame4 = mcstudy->plotNLL(Bins(40));
97
98 // Make some histograms from the parameter dataset
99 TH1 *hh_cor_a0_s1f = mcstudy->fitParDataSet().createHistogram("hh", a1, YVar(sig1frac));
100 TH1 *hh_cor_a0_a1 = mcstudy->fitParDataSet().createHistogram("hh", a0, YVar(a1));
101
102 // Access some of the saved fit results from individual toys
103 TH2 *corrHist000 = mcstudy->fitResult(0)->correlationHist("c000");
104 TH2 *corrHist127 = mcstudy->fitResult(127)->correlationHist("c127");
105 TH2 *corrHist953 = mcstudy->fitResult(953)->correlationHist("c953");
106
107 // Draw all plots on a canvas
108 gStyle->SetOptStat(0);
109 TCanvas *c = new TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900);
110 c->Divide(3, 3);
111 c->cd(1);
113 frame1->GetYaxis()->SetTitleOffset(1.4);
114 frame1->Draw();
115 c->cd(2);
117 frame2->GetYaxis()->SetTitleOffset(1.4);
118 frame2->Draw();
119 c->cd(3);
121 frame3->GetYaxis()->SetTitleOffset(1.4);
122 frame3->Draw();
123 c->cd(4);
125 frame4->GetYaxis()->SetTitleOffset(1.4);
126 frame4->Draw();
127 c->cd(5);
129 hh_cor_a0_s1f->GetYaxis()->SetTitleOffset(1.4);
130 hh_cor_a0_s1f->Draw("box");
131 c->cd(6);
133 hh_cor_a0_a1->GetYaxis()->SetTitleOffset(1.4);
134 hh_cor_a0_a1->Draw("box");
135 c->cd(7);
137 corrHist000->GetYaxis()->SetTitleOffset(1.4);
138 corrHist000->Draw("colz");
139 c->cd(8);
141 corrHist127->GetYaxis()->SetTitleOffset(1.4);
142 corrHist127->Draw("colz");
143 c->cd(9);
145 corrHist953->GetYaxis()->SetTitleOffset(1.4);
146 corrHist953->Draw("colz");
147
148 // Make RooMCStudy object available on command line after
149 // macro finishes
151}
#define c(i)
Definition RSha256.hxx:101
#define gDirectory
Definition TDirectory.h:384
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Efficient implementation of a sum of PDFs of the form.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Chebychev polynomial p.d.f.
TH2 * correlationHist(const char *name="correlation_matrix") const
Return TH2D of correlation matrix.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
Helper class to facilitate Monte Carlo studies such as 'goodness-of-fit' studies, that involve fittin...
Definition RooMCStudy.h:32
RooPlot * plotPull(const RooRealVar &param, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of pull values for the specified parameter on a newly created frame.
const RooFitResult * fitResult(Int_t sampleNum) const
Return the RooFitResult of the fit with the given run number.
RooPlot * plotNLL(const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the -log(L) values on a newly created frame.
const RooDataSet & fitParDataSet()
Return a RooDataSet containing the post-fit parameters of each toy cycle.
bool generateAndFit(Int_t nSamples, Int_t nEvtPerSample=0, bool keepGenData=false, const char *asciiFilePat=nullptr)
Generate and fit 'nSamples' samples of 'nEvtPerSample' events.
RooPlot * plotParam(const RooRealVar &param, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the fitted value of the given parameter on a newly created frame.
RooPlot * plotError(const RooRealVar &param, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={})
Plot the distribution of the fit errors for the specified parameter on a newly created frame.
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
TAxis * GetYaxis() const
Definition RooPlot.cxx:1276
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:649
Variable that can be changed from the outside.
Definition RooRealVar.h:37
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
The Canvas class.
Definition TCanvas.h:23
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TAxis * GetYaxis()
Definition TH1.h:325
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
Service class for 2-D histogram classes.
Definition TH2.h:30
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:1640
RooCmdArg FitGauss(bool flag=true)
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg={})
RooCmdArg Silence(bool flag=true)
RooCmdArg FitOptions(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={})
RooCmdArg Binned(bool flag=true)
RooCmdArg Bins(Int_t nbin)
RooCmdArg Save(bool flag=true)
RooCmdArg PrintEvalErrors(Int_t numErrors)
RooCmdArg Extended(bool flag=true)
Double_t x[n]
Definition legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26