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