ROOT   Reference Guide
rf613_global_observables.py
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) * 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) * 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) * Poisson(count_obs | count) \f]
48##
49## Unlike a Guassian, 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## \macro_code
54##
55## \date January 2022
56## \author Jonas Rembser
57
58
59import ROOT
60
61
62# Setting up the model and creating toy dataset
63# ---------------------------------------------
64
65# l'(x | mu, sigma) = l(x | mu, sigma) * Gauss(mu_obs | mu, 0.2)
66
67# event observables
68x = ROOT.RooRealVar("x", "x", -10, 10)
69
70# parameters
71mu = ROOT.RooRealVar("mu", "mu", 0.0, -10, 10)
72sigma = ROOT.RooRealVar("sigma", "sigma", 1.0, 0.1, 2.0)
73
74# Gaussian model for event observables
75gauss = ROOT.RooGaussian("gauss", "gauss", x, mu, sigma)
76
77# global observables (which are not parameters so they are constant)
78mu_obs = ROOT.RooRealVar("mu_obs", "mu_obs", 1.0, -10, 10)
79mu_obs.setConstant()
80# note: alternatively, one can create a constant with default limits using RooRealVar("mu_obs", "mu_obs", 1.0)
81
82# constraint pdf
83constraint = ROOT.RooGaussian("constraint", "constraint", mu_obs, mu, ROOT.RooFit.RooConst(0.2))
84
85# full pdf including constraint pdf
86model = ROOT.RooProdPdf("model", "model", [gauss, constraint])
87
88# Generating toy data with randomized global observables
89# ------------------------------------------------------
90
91# For most toy-based statistical procedures, it is necessary to also
92# randomize the global observable when generating toy datasets.
93#
94# To that end, let's generate a single event from the model and take the
95# global observable value (the same is done in the RooStats:ToyMCSampler
96# class):
97
98dataGlob = model.generate({mu_obs}, 1)
99
100# Next, we temporarily set the value of mu_obs to the randomized value for
101# generating our toy dataset:
102mu_obs_orig_val = mu_obs.getVal()
103
104ROOT.RooArgSet(mu_obs).assign(dataGlob.get(0))
105
106# actually generate the toy dataset
107data = model.generate({x}, 1000)
108
109# When fitting the toy dataset, it is important to set the global
110# observables in the fit to the values that were used to generate the toy
111# dataset. To facilitate the bookkeeping of global observable values, you
112# can attach a snapshot with the current global observable values to the
113# dataset like this (new feature introduced in ROOT 6.26):
114
115data.setGlobalObservables({mu_obs})
116
117# reset original mu_obs value
118mu_obs.setVal(mu_obs_orig_val)
119
120# Fitting a model with global observables
121# ---------------------------------------
122
123# Create snapshot of original parameters to reset parameters after fitting
124modelParameters = model.getParameters(data.get())
125origParameters = modelParameters.snapshot()
126
127# When you fit a model that includes global observables, you need to
128# specify them in the call to RooAbsPdf::fitTo with the
129# RooFit::GlobalObservables command argument. By default, the global
130# observable values attached to the dataset will be prioritized over the
131# values in the model, so the following fit correctly uses the randomized
132# global observable values from the toy dataset:
133print("1. model.fitTo(*data, GlobalObservables(mu_obs))")
134print("------------------------------------------------\n")
135model.fitTo(data, GlobalObservables=mu_obs, PrintLevel=-1, Save=True).Print()
136modelParameters.assign(origParameters)
137
138# In our example, the set of global observables is attached to the toy
139# dataset. In this case, you can actually drop the GlobalObservables()
140# command argument, because the global observables are automatically
141# figured out from the data set (this fit result should be identical to the
142# previous one).
143print("2. model.fitTo(*data)")
144print("---------------------\n")
145model.fitTo(data, PrintLevel=-1, Save=True).Print()
146modelParameters.assign(origParameters)
147
148# If you want to explicitly ignore the global observables in the dataset,
149# you can do that by specifying GlobalObservablesSource("model"). Keep in
150# mind that now it's also again your responsability to define the set of
151# global observables.
152print('3. model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource("model"))')
153print("------------------------------------------------\n")
154model.fitTo(data, GlobalObservables=mu_obs, GlobalObservablesSource="model", PrintLevel=-1, Save=True).Print()
155modelParameters.assign(origParameters)
void Print(std::ostream &os, const OptionType &opt)