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