Logo ROOT  
Reference Guide
rf709_BarlowBeeston.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 /// Implementing the Barlow-Beeston method for taking into account the statistical
5 /// uncertainty of a Monte-Carlo fit template.
6 ///
7 /// \macro_image
8 /// \macro_output
9 /// \macro_code
10 ///
11 /// Based on a demo by Wouter Verkerke
12 //
13 /// \date June 2019
14 /// \author Stephan Hageboeck, CERN
15 
16 #include "RooRealVar.h"
17 #include "RooGaussian.h"
18 #include "RooUniform.h"
19 #include "RooDataSet.h"
20 #include "RooDataHist.h"
21 #include "RooHistFunc.h"
22 #include "RooRealSumPdf.h"
23 #include "RooParamHistFunc.h"
24 #include "RooHistConstraint.h"
25 #include "RooProdPdf.h"
26 #include "RooPlot.h"
27 
28 #include "TCanvas.h"
29 #include "TPaveText.h"
30 
31 #include <iostream>
32 #include <memory>
33 
34 using namespace RooFit;
35 
36 void rf709_BarlowBeeston()
37 {
38  // First, construct a likelihood model with a Gaussian signal on top of a uniform background
39  RooRealVar x("x", "x", -20, 20);
40  x.setBins(25);
41 
42  RooRealVar meanG("meanG", "meanG", 1, -10, 10);
43  RooRealVar sigG("sigG", "sigG", 1.5, -10, 10);
44  RooGaussian g("g", "Gauss", x, meanG, sigG);
45  RooUniform u("u", "Uniform", x);
46 
47 
48  // Generate the data to be fitted
49  std::unique_ptr<RooDataSet> sigData(g.generate(x, 50));
50  std::unique_ptr<RooDataSet> bkgData(u.generate(x, 1000));
51 
52  RooDataSet sumData("sumData", "Gauss + Uniform", x);
53  sumData.append(*sigData);
54  sumData.append(*bkgData);
55 
56 
57  // Make histogram templates for signal and background.
58  // Let's take a signal distribution with low statistics and a more accurate
59  // background distribution.
60  // Normally, these come from Monte Carlo simulations, but we will just generate them.
61  std::unique_ptr<RooDataHist> dh_sig( g.generateBinned(x, 50) );
62  std::unique_ptr<RooDataHist> dh_bkg( u.generateBinned(x, 10000) );
63 
64 
65  // ***** Case 0 - 'Rigid templates' *****
66 
67  // Construct histogram shapes for signal and background
68  RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig);
69  RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg);
70 
71  // Construct scale factors for adding the two distributions
72  RooRealVar Asig0("Asig","Asig",1,0.01,5000);
73  RooRealVar Abkg0("Abkg","Abkg",1,0.01,5000);
74 
75  // Construct the sum model
76  RooRealSumPdf model0("model0","model0",
77  RooArgList(p_h_sig,p_h_bkg),
78  RooArgList(Asig0,Abkg0),
79  true);
80 
81 
82 
83  // ***** Case 1 - 'Barlow Beeston' *****
84 
85  // Construct parameterized histogram shapes for signal and background
86  RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig);
87  RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg);
88 
89  RooRealVar Asig1("Asig","Asig",1,0.01,5000);
90  RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000);
91 
92  // Construct the sum of these
93  RooRealSumPdf model_tmp("sp_ph", "sp_ph",
94  RooArgList(p_ph_sig1,p_ph_bkg1),
95  RooArgList(Asig1,Abkg1),
96  true);
97 
98  // Construct the subsidiary poisson measurements constraining the histogram parameters
99  // These ensure that the bin contents of the histograms are only allowed to vary within
100  // the statistical uncertainty of the Monte Carlo.
101  RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1);
102  RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1);
103 
104  // Construct the joint model with template PDFs and constraints
105  RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x));
106 
107 
108 
109  // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) *****
110 
111  // Construct the histogram shapes, using the same parameters for signal and background
112  // This requires passing the first histogram to the second, so that their common parameters
113  // can be re-used.
114  // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`.
115  // This allows bin 0 to fluctuate up and down.
116  // Then, the SAME parameters are connected to the background histogram, so the bins flucutate
117  // synchronously. This reduces the number of parameters.
118  RooParamHistFunc p_ph_sig2("p_ph_sig2", "p_ph_sig2", *dh_sig);
119  RooParamHistFunc p_ph_bkg2("p_ph_bkg2", "p_ph_bkg2", *dh_bkg, p_ph_sig2, true);
120 
121  RooRealVar Asig2("Asig","Asig",1,0.01,5000);
122  RooRealVar Abkg2("Abkg","Abkg",1,0.01,5000);
123 
124  // As before, construct the sum of signal2 and background2
125  RooRealSumPdf model2_tmp("sp_ph","sp_ph",
126  RooArgList(p_ph_sig2,p_ph_bkg2),
127  RooArgList(Asig2,Abkg2),
128  true);
129 
130  // Construct the subsidiary poisson measurements constraining the statistical fluctuations
131  RooHistConstraint hc_sigbkg("hc_sigbkg","hc_sigbkg",RooArgSet(p_ph_sig2,p_ph_bkg2));
132 
133  // Construct the joint model
134  RooProdPdf model2("model2","model2",hc_sigbkg, RooFit::Conditional(model2_tmp,x));
135 
136 
137 
138  // ************ Fit all models to data and plot *********************
139 
140  auto result0 = model0.fitTo(sumData, PrintLevel(0), Save());
141  auto result1 = model1.fitTo(sumData, PrintLevel(0), Save());
142  auto result2 = model2.fitTo(sumData, PrintLevel(0), Save());
143 
144 
145  TCanvas* can = new TCanvas("can", "", 1500, 600);
146  can->Divide(3,1);
147 
148  TPaveText pt(-19.5, 1, -2, 25);
149  pt.SetFillStyle(0);
150  pt.SetBorderSize(0);
151 
152 
153  can->cd(1);
154  auto frame = x.frame(Title("No template uncertainties"));
155  // Plot data to enable automatic determination of model0 normalisation:
156  sumData.plotOn(frame);
157  model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0));
158  // Plot data again to show it on top of model0 error bands:
159  sumData.plotOn(frame);
160  // Plot model components
161  model0.plotOn(frame, LineColor(kBlue));
162  model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure));
163  model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed));
164  model0.paramOn(frame);
165 
166  sigData->plotOn(frame, MarkerColor(kBlue));
167  frame->Draw();
168 
169  for (auto text : {
170  "No template uncertainties",
171  "are taken into account.",
172  "This leads to low errors",
173  "for the parameters A, since",
174  "the only source of errors",
175  "are the data statistics."}) {
176  pt.AddText(text);
177  }
178  pt.DrawClone();
179 
180 
181  can->cd(2);
182  frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately"));
183  sumData.plotOn(frame);
184  model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1));
185  // Plot data again to show it on top of error bands:
186  sumData.plotOn(frame);
187  model1.plotOn(frame, LineColor(kBlue));
188  model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure));
189  model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed));
190  model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1)));
191 
192  sigData->plotOn(frame, MarkerColor(kBlue));
193  frame->Draw();
194 
195  pt.Clear();
196  for (auto text : {
197  "With gamma parameters, the",
198  "signal & background templates",
199  "can adapt to the data.",
200  "Note how the blue signal",
201  "template changes its shape.",
202  "This leads to higher errors",
203  "of the scale parameters A."}) {
204  pt.AddText(text);
205  }
206  pt.DrawClone();
207 
208  can->cd(3);
209  frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)"));
210  sumData.plotOn(frame);
211  model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2));
212  // Plot data again to show it on top of model0 error bands:
213  sumData.plotOn(frame);
214  model2.plotOn(frame, LineColor(kBlue));
215  model2.plotOn(frame, Components(p_ph_sig2), LineColor(kAzure));
216  model2.plotOn(frame, Components(p_ph_bkg2), LineColor(kRed));
217  model2.paramOn(frame, Parameters(RooArgSet(Asig2, Abkg2)));
218 
219  sigData->plotOn(frame, MarkerColor(kBlue));
220  frame->Draw();
221 
222  pt.Clear();
223  for (auto text : {
224  "When signal and background",
225  "template share one gamma para-",
226  "meter per bin, they adapt less.",
227  "The errors of the A parameters",
228  "also shrink slightly."}) {
229  pt.AddText(text);
230  }
231  pt.DrawClone();
232 
233 
234  std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl;
235  std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl;
236  std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl;
237 
238 }
TPave::SetBorderSize
virtual void SetBorderSize(Int_t bordersize=4)
Definition: TPave.h:73
RooUniform
Flat p.d.f.
Definition: RooUniform.h:24
TPaveText::Clear
virtual void Clear(Option_t *option="")
Clear all lines in this pavetext.
Definition: TPaveText.cxx:209
TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
RooHistFunc.h
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition: RooGlobalFunc.cxx:192
RooFit::VisualizeError
RooCmdArg VisualizeError(const RooDataSet &paramData, Double_t Z=1)
Definition: RooGlobalFunc.cxx:71
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooGaussian.h
TObject::DrawClone
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad for instance with: gROOT->SetSelectedPad(gPad...
Definition: TObject.cxx:221
x
Double_t x[n]
Definition: legend1.C:17
RooGaussian
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
TPad::Divide
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1177
TCanvas.h
RooDataSet.h
text
TText * text
Definition: entrylist_figure1.C:10
TPaveText.h
RooProdPdf.h
RooParamHistFunc.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:708
RooDataHist.h
RooPlot.h
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooHistConstraint
The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
Definition: RooHistConstraint.h:19
RooHistFunc
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
Definition: RooHistFunc.h:30
TPaveText::AddText
virtual TText * AddText(Double_t x1, Double_t y1, const char *label)
Add a new Text line to this pavetext at given coordinates.
Definition: TPaveText.cxx:183
RooFit::Conditional
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, Bool_t depsAreCond=kFALSE)
Definition: RooGlobalFunc.cxx:231
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:57
RooHistConstraint.h
TCanvas
The Canvas class.
Definition: TCanvas.h:23
RooParamHistFunc
A histogram function that assigns scale parameters to every bin.
Definition: RooParamHistFunc.h:20
RooFit::Parameters
RooCmdArg Parameters(const RooArgSet &params)
Definition: RooGlobalFunc.cxx:224
kBlue
@ kBlue
Definition: Rtypes.h:66
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
TPaveText
A Pave (see TPave) with text, lines or/and boxes inside.
Definition: TPaveText.h:21
RooFit::MarkerColor
RooCmdArg MarkerColor(Color_t color)
Definition: RooGlobalFunc.cxx:88
pt
TPaveText * pt
Definition: entrylist_figure1.C:7
RooRealSumPdf.h
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:37
RooFit::Components
RooCmdArg Components(const RooArgSet &compSet)
Definition: RooGlobalFunc.cxx:75
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:190
RooProdPdf
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:37
kAzure
@ kAzure
Definition: Rtypes.h:67
RooRealSumPdf
The class RooRealSumPdf implements a PDF constructed from a sum of functions:
Definition: RooRealSumPdf.h:25
RooFit::Title
RooCmdArg Title(const char *name)
Definition: RooGlobalFunc.cxx:176
RooUniform.h
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
g
#define g(i)
Definition: RSha256.hxx:105