Logo ROOT  
Reference Guide
rf611_weightedfits.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 /// Likelihood and minimization: Parameter uncertainties for weighted unbinned ML fits
5 ///
6 /// ## Parameter uncertainties for weighted unbinned ML fits
7 ///
8 /// Based on example from https://arxiv.org/abs/1911.01303
9 ///
10 /// This example compares different approaches to determining parameter uncertainties in weighted unbinned maximum likelihood fits.
11 /// Performing a weighted unbinned maximum likelihood fits can be useful to account for acceptance effects and to statistically subtract background events using the sPlot formalism.
12 /// It is however well known that the inverse Hessian matrix does not yield parameter uncertainties with correct coverage in the presence of event weights.
13 /// Three approaches to the determination of parameter uncertainties are compared in this example:
14 ///
15 /// 1. Using the inverse weighted Hessian matrix [`SumW2Error(false)`]
16 ///
17 /// 2. Using the expression [`SumW2Error(true)`]
18 /// \f[
19 /// V_{ij} = H_{ik}^{-1} C_{kl} H_{lj}^{-1}
20 /// \f]
21 /// where H is the weighted Hessian matrix and C is the Hessian matrix with squared weights
22 ///
23 /// 3. The asymptotically correct approach (for details please see https://arxiv.org/abs/1911.01303) [`Asymptotic(true)`]
24 /// \f[
25 /// V_{ij} = H_{ik}^{-1} D_{kl} H_{lj}^{-1}
26 /// \f]
27 /// where H is the weighted Hessian matrix and D is given by
28 /// \f[
29 /// D_{kl} = \sum_{e=1}^{N} w_e^2 \frac{\partial \log(P)}{\partial \lambda_k}\frac{\partial \log(P)}{\partial \lambda_l}
30 /// \f]
31 /// with the event weight \f$w_e\f$.
32 ///
33 /// The example performs the fit of a second order polynomial in the angle cos(theta) [-1,1] to a weighted data set.
34 /// The polynomial is given by
35 /// \f[
36 /// P = \frac{ 1 + c_0 \cdot \cos(\theta) + c_1 \cdot \cos(\theta) \cdot \cos(\theta) }{\mathrm{Norm}}
37 /// \f]
38 /// The two coefficients \f$ c_0 \f$ and \f$ c_1 \f$ and their uncertainties are to be determined in the fit.
39 ///
40 /// The per-event weight is used to correct for an acceptance effect, two different acceptance models can be studied:
41 /// - `acceptancemodel==1`: eff = \f$ 0.3 + 0.7 \cdot \cos(\theta) \cdot \cos(\theta) \f$
42 /// - `acceptancemodel==2`: eff = \f$ 1.0 - 0.7 \cdot \cos(\theta) \cdot \cos(\theta) \f$
43 /// The data is generated to be flat before the acceptance effect.
44 ///
45 /// The performance of the different approaches to determine parameter uncertainties is compared using the pull distributions from a large number of pseudoexperiments.
46 /// The pull is defined as \f$ (\lambda_i - \lambda_{gen})/\sigma(\lambda_i) \f$, where \f$ \lambda_i \f$ is the fitted parameter and \f$ \sigma(\lambda_i) \f$ its
47 /// uncertainty for pseudoexperiment number i.
48 /// If the fit is unbiased and the parameter uncertainties are estimated correctly, the pull distribution should be a Gaussian centered around zero with a width of one.
49 ///
50 /// \macro_image
51 /// \macro_output
52 /// \macro_code
53 ///
54 /// \date November 2019
55 /// \author Christoph Langenbruch
56 
57 #include "TH1D.h"
58 #include "TCanvas.h"
59 #include "TROOT.h"
60 #include "TStyle.h"
61 #include "TRandom3.h"
62 #include "TLegend.h"
63 #include "RooRealVar.h"
64 #include "RooFitResult.h"
65 #include "RooDataSet.h"
66 #include "RooPolynomial.h"
67 
68 using namespace RooFit;
69 
70 
71 int rf611_weightedfits(int acceptancemodel=2) {
72  // I n i t i a l i s a t i o n a n d S e t u p
73  //------------------------------------------------
74 
75  //plotting options
76  gStyle->SetPaintTextFormat(".1f");
77  gStyle->SetEndErrorSize(6.0);
78  gStyle->SetTitleSize(0.05, "XY");
79  gStyle->SetLabelSize(0.05, "XY");
80  gStyle->SetTitleOffset(0.9, "XY");
81  gStyle->SetTextSize(0.05);
82  gStyle->SetPadLeftMargin(0.125);
83  gStyle->SetPadBottomMargin(0.125);
84  gStyle->SetPadTopMargin(0.075);
85  gStyle->SetPadRightMargin(0.075);
87  gStyle->SetMarkerSize(1.0);
90 
91  //initialise TRandom3
92  TRandom3* rnd = new TRandom3();
93  rnd->SetSeed(191101303);
94 
95  //accepted events and events weighted to account for the acceptance
96  TH1D* haccepted = new TH1D("haccepted", "Generated events;cos(#theta);#events", 40, -1.0, 1.0);
97  TH1D* hweighted = new TH1D("hweighted", "Generated events;cos(#theta);#events", 40, -1.0, 1.0);
98  //histograms holding pull distributions
99  //using the inverse Hessian matrix
100  TH1D* hc0pull1 = new TH1D("hc0pull1", "Inverse weighted Hessian matrix [SumW2Error(false)];Pull (c_{0}^{fit}-c_{0}^{gen})/#sigma(c_{0});", 20, -5.0, 5.0);
101  TH1D* hc1pull1 = new TH1D("hc1pull1", "Inverse weighted Hessian matrix [SumW2Error(false)];Pull (c_{1}^{fit}-c_{1}^{gen})/#sigma(c_{1});", 20, -5.0, 5.0);
102  //using the correction with the Hessian matrix with squared weights
103  TH1D* hc0pull2 = new TH1D("hc0pull2", "Hessian matrix with squared weights [SumW2Error(true)];Pull (c_{0}^{fit}-c_{0}^{gen})/#sigma(c_{0});", 20, -5.0, 5.0);
104  TH1D* hc1pull2 = new TH1D("hc1pull2", "Hessian matrix with squared weights [SumW2Error(true)];Pull (c_{1}^{fit}-c_{1}^{gen})/#sigma(c_{1});", 20, -5.0, 5.0);
105  //asymptotically correct approach
106  TH1D* hc0pull3 = new TH1D("hc0pull3", "Asymptotically correct approach [Asymptotic(true)];Pull (c_{0}^{fit}-c_{0}^{gen})/#sigma(c_{0});", 20, -5.0, 5.0);
107  TH1D* hc1pull3 = new TH1D("hc1pull3", "Asymptotically correct approach [Asymptotic(true)];Pull (c_{1}^{fit}-c_{1}^{gen})/#sigma(c_{1});", 20, -5.0, 5.0);
108 
109  //number of pseudoexperiments (toys) and number of events per pseudoexperiment
110  constexpr unsigned int ntoys = 500;
111  constexpr unsigned int nstats = 5000;
112  //parameters used in the generation
113  constexpr double c0gen = 0.0;
114  constexpr double c1gen = 0.0;
115 
116  // Silence fitting and minimisation messages
117  auto& msgSv = RooMsgService::instance();
118  msgSv.getStream(1).removeTopic(RooFit::Minimization);
119  msgSv.getStream(1).removeTopic(RooFit::Fitting);
120 
121  std::cout << "Running " << ntoys*3 << " toy fits ..." << std::endl;
122 
123  // M a i n l o o p : r u n p s e u d o e x p e r i m e n t s
124  //----------------------------------------------------------------
125  for (unsigned int i=0; i<ntoys; i++) {
126  //S e t u p p a r a m e t e r s a n d P D F
127  //-----------------------------------------------
128  //angle theta and the weight to account for the acceptance effect
129  RooRealVar costheta("costheta","costheta", -1.0, 1.0);
130  RooRealVar weight("weight","weight", 0.0, 1000.0);
131 
132  //initialise parameters to fit
133  RooRealVar c0("c0","0th-order coefficient", c0gen, -1.0, 1.0);
134  RooRealVar c1("c1","1st-order coefficient", c1gen, -1.0, 1.0);
135  c0.setError(0.01);
136  c1.setError(0.01);
137  //create simple second-order polynomial as probability density function
138  RooPolynomial pol("pol", "pol", costheta, RooArgList(c0, c1), 1);
139 
140  //G e n e r a t e d a t a s e t f o r p s e u d o e x p e r i m e n t i
141  //-------------------------------------------------------------------------------
142  RooDataSet data("data","data",RooArgSet(costheta, weight), WeightVar("weight"));
143  //generate nstats events
144  for (unsigned int j=0; j<nstats; j++) {
145  bool finished = false;
146  //use simple accept/reject for generation
147  while (!finished) {
148  costheta = 2.0*rnd->Rndm()-1.0;
149  //efficiency for the specific value of cos(theta)
150  double eff = 1.0;
151  if (acceptancemodel == 1)
152  eff = 1.0 - 0.7 * costheta.getVal()*costheta.getVal();
153  else
154  eff = 0.3 + 0.7 * costheta.getVal()*costheta.getVal();
155  //use 1/eff as weight to account for acceptance
156  weight = 1.0/eff;
157  //accept/reject
158  if (10.0*rnd->Rndm() < eff*pol.getVal())
159  finished = true;
160  }
161  haccepted->Fill(costheta.getVal());
162  hweighted->Fill(costheta.getVal(), weight.getVal());
163  data.add(RooArgSet(costheta, weight), weight.getVal());
164  }
165 
166  //F i t t o y u s i n g t h e t h r e e d i f f e r e n t a p p r o a c h e s t o u n c e r t a i n t y d e t e r m i n a t i o n
167  //-------------------------------------------------------------------------------------------------------------------------------------------------
168  //this uses the inverse weighted Hessian matrix
169  RooFitResult* result = pol.fitTo(data, Save(true), SumW2Error(false), PrintLevel(-1), BatchMode(true));
170  hc0pull1->Fill((c0.getVal()-c0gen)/c0.getError());
171  hc1pull1->Fill((c1.getVal()-c1gen)/c1.getError());
172 
173  //this uses the correction with the Hesse matrix with squared weights
174  result = pol.fitTo(data, Save(true), SumW2Error(true), PrintLevel(-1), BatchMode(true));
175  hc0pull2->Fill((c0.getVal()-c0gen)/c0.getError());
176  hc1pull2->Fill((c1.getVal()-c1gen)/c1.getError());
177 
178  //this uses the asymptotically correct approach
179  result = pol.fitTo(data, Save(true), AsymptoticError(true), PrintLevel(-1), BatchMode(true));
180  hc0pull3->Fill((c0.getVal()-c0gen)/c0.getError());
181  hc1pull3->Fill((c1.getVal()-c1gen)/c1.getError());
182  }
183 
184  std::cout << "... done." << std::endl;
185 
186  // P l o t o u t p u t d i s t r i b u t i o n s
187  //--------------------------------------------------
188 
189  //plot accepted (weighted) events
190  gStyle->SetOptStat(0);
191  gStyle->SetOptFit(0);
192  TCanvas* cevents = new TCanvas("cevents", "cevents", 800, 600);
193  cevents->cd(1);
194  hweighted->SetMinimum(0.0);
195  hweighted->SetLineColor(2);
196  hweighted->Draw("hist");
197  haccepted->Draw("same hist");
198  TLegend* leg = new TLegend(0.6, 0.8, 0.9, 0.9);
199  leg->AddEntry(haccepted, "Accepted");
200  leg->AddEntry(hweighted, "Weighted");
201  leg->Draw();
202  cevents->Update();
203 
204  //plot pull distributions
205  TCanvas* cpull = new TCanvas("cpull", "cpull", 1200, 800);
206  cpull->Divide(3,2);
207  cpull->cd(1);
208  gStyle->SetOptStat(1100);
209  gStyle->SetOptFit(11);
210  hc0pull1->Fit("gaus");
211  hc0pull1->Draw("ep");
212  cpull->cd(2);
213  hc0pull2->Fit("gaus");
214  hc0pull2->Draw("ep");
215  cpull->cd(3);
216  hc0pull3->Fit("gaus");
217  hc0pull3->Draw("ep");
218  cpull->cd(4);
219  hc1pull1->Fit("gaus");
220  hc1pull1->Draw("ep");
221  cpull->cd(5);
222  hc1pull2->Fit("gaus");
223  hc1pull2->Draw("ep");
224  cpull->cd(6);
225  hc1pull3->Fit("gaus");
226  hc1pull3->Draw("ep");
227  cpull->Update();
228 
229  return 0;
230 }
RooFit::BatchMode
RooCmdArg BatchMode(bool flag=true)
Definition: RooGlobalFunc.cxx:159
TStyle::SetLabelSize
void SetLabelSize(Float_t size=0.04, Option_t *axis="X")
Set size of axis labels.
Definition: TStyle.cxx:1392
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
RooFit::Minimization
@ Minimization
Definition: RooGlobalFunc.h:67
TStyle::SetPadBottomMargin
void SetPadBottomMargin(Float_t margin=0.1)
Definition: TStyle.h:341
TStyle::SetPaintTextFormat
void SetPaintTextFormat(const char *format="g")
Definition: TStyle.h:369
TH1::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:399
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition: RooGlobalFunc.cxx:192
TStyle::SetTitleSize
void SetTitleSize(Float_t size=0.02, Option_t *axis="X")
Definition: TStyle.cxx:1767
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:210
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
TLegend.h
RooFit::WeightVar
RooCmdArg WeightVar(const char *name, Bool_t reinterpretAsWeight=kFALSE)
Definition: RooGlobalFunc.cxx:124
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
TStyle.h
TStyle::SetEndErrorSize
void SetEndErrorSize(Float_t np=2)
Set the size (in pixels) of the small lines drawn at the end of the error bars (TH1 or TGraphErrors).
Definition: TStyle.cxx:1288
TAttText::SetTextSize
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition: TAttText.h:46
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
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TCanvas.h
RooDataSet.h
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
TStyle::SetOptFit
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1541
RooPolynomial.h
TROOT.h
gStyle
R__EXTERN TStyle * gStyle
Definition: TStyle.h:412
TRandom3
Random number generator class based on M.
Definition: TRandom3.h:27
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
TRandom3::SetSeed
virtual void SetSeed(ULong_t seed=0)
Set the random generator sequence if seed is 0 (default value) a TUUID is generated and used to fill ...
Definition: TRandom3.cxx:206
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:708
RooPolynomial
RooPolynomial implements a polynomial p.d.f of the form.
Definition: RooPolynomial.h:28
RooFit::Fitting
@ Fitting
Definition: RooGlobalFunc.h:67
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3327
TRandom3.h
TStyle::SetHistLineWidth
void SetHistLineWidth(Width_t width=1)
Definition: TStyle.h:366
TStyle::SetPadTopMargin
void SetPadTopMargin(Float_t margin=0.1)
Definition: TStyle.h:342
TStyle::SetPadLeftMargin
void SetPadLeftMargin(Float_t margin=0.1)
Definition: TStyle.h:343
RooRealVar.h
TStyle::SetOptStat
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:1589
RooFitResult.h
TStyle::SetPadRightMargin
void SetPadRightMargin(Float_t margin=0.1)
Definition: TStyle.h:344
TStyle::SetHistLineColor
void SetHistLineColor(Color_t color=1)
Definition: TStyle.h:363
TCanvas
The Canvas class.
Definition: TCanvas.h:23
TRandom3::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:99
TCanvas::Update
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2504
leg
leg
Definition: legend1.C:34
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
TH1::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3869
TLegend
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:37
TH1D.h
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:190
TStyle::SetTitleOffset
void SetTitleOffset(Float_t offset=1, Option_t *axis="X")
Specify a parameter offset to control the distance between the axis and the axis title.
Definition: TStyle.cxx:1748
RooMsgService::instance
static RooMsgService & instance()
Return reference to singleton instance.
Definition: RooMsgService.cxx:366
RooFit::AsymptoticError
RooCmdArg AsymptoticError(Bool_t flag)
Definition: RooGlobalFunc.cxx:211
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3050
c1
return c1
Definition: legend1.C:41