Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf613_global_observables.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4/// This tutorial explains the concept of global observables in RooFit, and
5/// showcases how their values can be stored either in the model or in the
6/// dataset.
7///
8/// ### Introduction
9///
10/// Note: in this tutorial, we are multiplying the likelihood with an additional
11/// likelihood to constrain the parameters with auxiliary measurements. This is
12/// different from the `rf604_constraints` tutorial, where the likelihood is
13/// multiplied with a Bayesian prior to constrain the parameters.
14///
15///
16/// With RooFit, you usually optimize some model parameters `p` to maximize the
17/// likelihood `L` given the per-event or per-bin observations `x`:
18///
19/// \f[ L( x | p ) \f]
20///
21/// Often, the parameters are constrained with some prior likelihood `C`, which
22/// doesn't depend on the observables `x`:
23///
24/// \f[ L'( x | p ) = L( x | p ) * C( p ) \f]
25///
26/// Usually, these constraint terms depend on some auxiliary measurements of
27/// other observables `g`. The constraint term is then the likelihood of the
28/// so-called global observables:
29///
30/// \f[ L'( x | p ) = L( x | p ) * C( g | p ) \f]
31///
32/// For example, think of a model where the true luminosity `lumi` is a
33/// nuisance parameter that is constrained by an auxiliary measurement
34/// `lumi_obs` with uncertainty `lumi_obs_sigma`:
35///
36/// \f[ L'(data | mu, lumi) = L(data | mu, lumi) * \text{Gauss}(lumi_obs | lumi, lumi_obs_sigma) \f]
37///
38/// As a Gaussian is symmetric under exchange of the observable and the mean
39/// parameter, you can also sometimes find this equivalent but less conventional
40/// formulation for Gaussian constraints:
41///
42/// \f[ L'(data | mu, lumi) = L(data | mu, lumi) * \text{Gauss}(lumi | lumi_obs, lumi_obs_sigma) \f]
43///
44/// If you wanted to constrain a parameter that represents event counts, you
45/// would use a Poissonian constraint, e.g.:
46///
47/// \f[ L'(data | mu, count) = L(data | mu, count) * \text{Poisson}(count_obs | count) \f]
48///
49/// Unlike a Gaussian, a Poissonian is not symmetric under exchange of the
50/// observable and the parameter, so here you need to be more careful to follow
51/// the global observable prescription correctly.
52///
53///
54/// \macro_code
55/// \macro_output
56///
57/// \date January 2022
58/// \author Jonas Rembser
59
60#include <RooArgSet.h>
61#include <RooConstVar.h>
62#include <RooDataSet.h>
63#include <RooFitResult.h>
64#include <RooGaussian.h>
65#include <RooProdPdf.h>
66#include <RooRealVar.h>
67
69{
70 using namespace RooFit;
71
72 // Silence info output for this tutorial
75
76 // Setting up the model and creating toy dataset
77 // ---------------------------------------------
78
79 // l'(x | mu, sigma) = l(x | mu, sigma) * Gauss(mu_obs | mu, 0.2)
80
81 // event observables
82 RooRealVar x("x", "x", -10, 10);
83
84 // parameters
85 RooRealVar mu("mu", "mu", 0.0, -10, 10);
86 RooRealVar sigma("sigma", "sigma", 1.0, 0.1, 2.0);
87
88 // Gaussian model for event observables
89 RooGaussian gauss("gauss", "gauss", x, mu, sigma);
90
91 // global observables (which are not parameters so they are constant)
92 RooRealVar mu_obs("mu_obs", "mu_obs", 1.0, -10, 10);
93 mu_obs.setConstant();
94 // note: alternatively, one can create a constant with default limits using `RooRealVar("mu_obs", "mu_obs", 1.0)`
95
96 // constraint pdf
97 RooGaussian constraint("constraint", "constraint", mu_obs, mu, 0.1);
98
99 // full pdf including constraint pdf
100 RooProdPdf model("model", "model", {gauss, constraint});
101
102 // Generating toy data with randomized global observables
103 // ------------------------------------------------------
104
105 // For most toy-based statistical procedures, it is necessary to also
106 // randomize the global observable when generating toy datasets.
107 //
108 // To that end, let's generate a single event from the model and take the
109 // global observable value (the same is done in the RooStats:ToyMCSampler
110 // class):
111
112 std::unique_ptr<RooDataSet> dataGlob{model.generate({mu_obs}, 1)};
113
114 // Next, we temporarily set the value of `mu_obs` to the randomized value for
115 // generating our toy dataset:
116 double mu_obs_orig_val = mu_obs.getVal();
117
118 RooArgSet{mu_obs}.assign(*dataGlob->get(0));
119
120 // Actually generate the toy dataset. We don't generate too many events,
121 // otherwise, the constraint will not have much weight in the fit and the
122 // result looks like it's unaffected by it.
123 std::unique_ptr<RooDataSet> data{model.generate({x}, 50)};
124
125 // When fitting the toy dataset, it is important to set the global
126 // observables in the fit to the values that were used to generate the toy
127 // dataset. To facilitate the bookkeeping of global observable values, you
128 // can attach a snapshot with the current global observable values to the
129 // dataset like this (new feature introduced in ROOT 6.26):
130
131 data->setGlobalObservables({mu_obs});
132
133 // reset original mu_obs value
134 mu_obs.setVal(mu_obs_orig_val);
135
136 // Fitting a model with global observables
137 // ---------------------------------------
138
139 // Create snapshot of original parameters to reset parameters after fitting
140 RooArgSet modelParameters;
141 model.getParameters(data->get(), modelParameters);
142 RooArgSet origParameters;
143 modelParameters.snapshot(origParameters);
144
145 using FitRes = std::unique_ptr<RooFitResult>;
146
147 // When you fit a model that includes global observables, you need to
148 // specify them in the call to RooAbsPdf::fitTo with the
149 // RooFit::GlobalObservables command argument. By default, the global
150 // observable values attached to the dataset will be prioritized over the
151 // values in the model, so the following fit correctly uses the randomized
152 // global observable values from the toy dataset:
153 std::cout << "1. model.fitTo(*data, GlobalObservables(mu_obs))\n";
154 std::cout << "------------------------------------------------\n";
155 FitRes res1{model.fitTo(*data, GlobalObservables(mu_obs), PrintLevel(-1), Save())};
156 res1->Print();
157 modelParameters.assign(origParameters);
158
159 // In our example, the set of global observables is attached to the toy
160 // dataset. In this case, you can actually drop the GlobalObservables()
161 // command argument, because the global observables are automatically
162 // figured out from the data set (this fit result should be identical to the
163 // previous one).
164 std::cout << "2. model.fitTo(*data)\n";
165 std::cout << "---------------------\n";
166 FitRes res2{model.fitTo(*data, PrintLevel(-1), Save())};
167 res2->Print();
168 modelParameters.assign(origParameters);
169
170 // If you want to explicitly ignore the global observables in the dataset,
171 // you can do that by specifying GlobalObservablesSource("model"). Keep in
172 // mind that now it's also again your responsibility to define the set of
173 // global observables.
174 std::cout << "3. model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource(\"model\"))\n";
175 std::cout << "------------------------------------------------\n";
176 FitRes res3{model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource("model"), PrintLevel(-1), Save())};
177 res3->Print();
178 modelParameters.assign(origParameters);
179}
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:57
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:154
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
static RooMsgService & instance()
Return reference to singleton instance.
StreamConfig & getStream(Int_t id)
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
Variable that can be changed from the outside.
Definition RooRealVar.h:37
RooCmdArg Save(bool flag=true)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg GlobalObservablesSource(const char *sourceName)
RooCmdArg PrintLevel(Int_t code)
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
void removeTopic(RooFit::MsgTopic oldTopic)