Logo ROOT  
Reference Guide
rf802_mcstudy_addons.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4/// Validation and MC studies:
5/// RooMCStudy - using separate fit and generator models, using the chi^2 calculator model
6/// Running a biased fit model against an optimal fit.
7///
8/// \macro_image
9/// \macro_output
10/// \macro_code
11///
12/// \date July 2008
13/// \author Wouter Verkerke
14
15#include "RooRealVar.h"
16#include "RooDataSet.h"
17#include "RooGaussian.h"
18#include "RooConstVar.h"
19#include "RooChebychev.h"
20#include "RooAddPdf.h"
21#include "RooMCStudy.h"
22#include "RooChi2MCSModule.h"
23#include "RooPlot.h"
24#include "TCanvas.h"
25#include "TAxis.h"
26#include "TH1.h"
27#include "TDirectory.h"
28#include "TLegend.h"
29
30using namespace RooFit;
31
32void rf802_mcstudy_addons()
33{
34
35 // C r e a t e m o d e l
36 // -----------------------
37
38 // Observables, parameters
39 RooRealVar x("x", "x", -10, 10);
40 x.setBins(10);
41 RooRealVar mean("mean", "mean of gaussian", 0, -2., 1.8);
42 RooRealVar sigma("sigma", "width of gaussian", 5, 1, 10);
43
44 // Create Gaussian pdf
45 RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
46
47 // C r e a t e m a n a g e r w i t h c h i ^ 2 a d d - o n m o d u l e
48 // ----------------------------------------------------------------------------
49
50 // Create study manager for binned likelihood fits of a Gaussian pdf in 10 bins
51 RooMCStudy *mcs = new RooMCStudy(gauss, x, Silence(), Binned());
52
53 // Add chi^2 calculator module to mcs
54 RooChi2MCSModule chi2mod;
55 mcs->addModule(chi2mod);
56
57 // Generate 1000 samples of 1000 events
58 mcs->generateAndFit(2000, 1000);
59
60 // Fill histograms with distributions chi2 and prob(chi2,ndf) that
61 // are calculated by RooChiMCSModule
62 TH1 *hist_chi2 = mcs->fitParDataSet().createHistogram("chi2");
63 hist_chi2->SetTitle("#chi^{2} values of all toy runs;#chi^{2}");
64 TH1 *hist_prob = mcs->fitParDataSet().createHistogram("prob");
65 hist_prob->SetTitle("Corresponding #chi^{2} probability;Prob(#chi^{2},ndof)");
66
67
68 // C r e a t e m a n a g e r w i t h s e p a r a t e f i t m o d e l
69 // ----------------------------------------------------------------------------
70
71 // Create alternate pdf with shifted mean
72 RooRealVar mean2("mean2", "mean of gaussian 2", 2.);
73 RooGaussian gauss2("gauss2", "gaussian PDF2", x, mean2, sigma);
74
75 // Create study manager with separate generation and fit model. This configuration
76 // is set up to generate biased fits as the fit and generator model have different means,
77 // and the mean parameter is limited to [-2., 1.8], so it just misses the optimal
78 // mean value of 2 in the data.
79 RooMCStudy *mcs2 = new RooMCStudy(gauss2, x, FitModel(gauss), Silence(), Binned());
80
81 // Add chi^2 calculator module to mcs
82 RooChi2MCSModule chi2mod2;
83 mcs2->addModule(chi2mod2);
84
85 // Generate 1000 samples of 1000 events
86 mcs2->generateAndFit(2000, 1000);
87
88 // Request a the pull plot of mean. The pulls will be one-sided because
89 // `mean` is limited to 1.8.
90 // Note that RooFit will have trouble to compute the pulls because the parameters
91 // are called `mean` in the fit, but `mean2` in the generator model. It is not obvious
92 // that these are related. RooFit will nevertheless compute pulls, but complain that
93 // this is risky.
94 auto pullMeanFrame = mcs2->plotPull(mean);
95
96 // Fill histograms with distributions chi2 and prob(chi2,ndf) that
97 // are calculated by RooChiMCSModule
98 TH1 *hist2_chi2 = mcs2->fitParDataSet().createHistogram("chi2");
99 TH1 *hist2_prob = mcs2->fitParDataSet().createHistogram("prob");
100 hist2_chi2->SetLineColor(kRed);
101 hist2_prob->SetLineColor(kRed);
102
103 TLegend leg;
104 leg.AddEntry(hist_chi2, "Optimal fit", "L");
105 leg.AddEntry(hist2_chi2, "Biased fit", "L");
106 leg.SetBorderSize(0);
107 leg.SetFillStyle(0);
108
109 TCanvas *c = new TCanvas("rf802_mcstudy_addons", "rf802_mcstudy_addons", 800, 400);
110 c->Divide(3);
111 c->cd(1);
112 gPad->SetLeftMargin(0.15);
113 hist_chi2->GetYaxis()->SetTitleOffset(1.4);
114 hist_chi2->Draw();
115 hist2_chi2->Draw("esame");
116 leg.DrawClone();
117 c->cd(2);
118 gPad->SetLeftMargin(0.15);
119 hist_prob->GetYaxis()->SetTitleOffset(1.4);
120 hist_prob->Draw();
121 hist2_prob->Draw("esame");
122 c->cd(3);
123 pullMeanFrame->Draw();
124
125
126 // Make RooMCStudy object available on command line after
127 // macro finishes
128 gDirectory->Add(mcs);
129}
#define c(i)
Definition: RSha256.hxx:101
@ kRed
Definition: Rtypes.h:66
#define gDirectory
Definition: TDirectory.h:290
#define gPad
Definition: TVirtualPad.h:287
RooChi2MCSModule is an add-on module to RooMCStudy that calculates the chi-squared of fitted p....
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.
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
RooMCStudy is a helper class to facilitate Monte Carlo studies such as 'goodness-of-fit' studies,...
Definition: RooMCStudy.h:32
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.
const RooDataSet & fitParDataSet()
Return a RooDataSet containing the post-fit parameters of each toy cycle.
Definition: RooMCStudy.cxx:981
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
void addModule(RooAbsMCStudyModule &module)
Insert given RooMCStudy add-on module to the processing chain of this MCStudy object.
Definition: RooMCStudy.cxx:445
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:293
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
The Canvas class.
Definition: TCanvas.h:23
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6678
TAxis * GetYaxis()
Definition: TH1.h:321
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
RooCmdArg Binned(Bool_t flag=kTRUE)
RooCmdArg Silence(Bool_t flag=kTRUE)
const Double_t sigma
Double_t x[n]
Definition: legend1.C:17
leg
Definition: legend1.C:34
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
void FitModel(RooWorkspace *, std::string data_name="obsData")
static constexpr double gauss