ROOT   Reference Guide
rf613_global_observables.C File Reference

## Detailed Description

This tutorial explains the concept of global observables in RooFit, and showcases how their values can be stored either in the model or in the dataset.

### Introduction

Note: in this tutorial, we are multiplying the likelihood with an additional likelihood to constrain the parameters with auxiliary measurements. This is different from the rf604_constraints tutorial, where the likelihood is multiplied with a Bayesian prior to constrain the parameters.

With RooFit, you usually optimize some model parameters p to maximize the likelihood L given the per-event or per-bin observations x:

$L( x | p )$

Often, the parameters are constrained with some prior likelihood C, which doesn't depend on the observables x:

$L'( x | p ) = L( x | p ) * C( p )$

Usually, these constraint terms depend on some auxiliary measurements of other observables g. The constraint term is then the likelihood of the so-called global observables:

$L'( x | p ) = L( x | p ) * C( g | p )$

For example, think of a model where the true luminosity lumi is a nuisance parameter that is constrained by an auxiliary measurement lumi_obs with uncertainty lumi_obs_sigma:

$L'(data | mu, lumi) = L(data | mu, lumi) * \text{Gauss}(lumi_obs | lumi, lumi_obs_sigma)$

As a Gaussian is symmetric under exchange of the observable and the mean parameter, you can also sometimes find this equivalent but less conventional formulation for Gaussian constraints:

$L'(data | mu, lumi) = L(data | mu, lumi) * \text{Gauss}(lumi | lumi_obs, lumi_obs_sigma)$

If you wanted to constrain a parameter that represents event counts, you would use a Poissonian constraint, e.g.:

$L'(data | mu, count) = L(data | mu, count) * \text{Poisson}(count_obs | count)$

Unlike a Guassian, a Poissonian is not symmetric under exchange of the observable and the parameter, so here you need to be more careful to follow the global observable prescription correctly.

#include <RooArgSet.h>
#include <RooConstVar.h>
#include <RooDataSet.h>
#include <RooFitResult.h>
#include <RooGaussian.h>
#include <RooProdPdf.h>
#include <RooRealVar.h>
{
using namespace RooFit;
// Silence info output for this tutorial
// Setting up the model and creating toy dataset
// ---------------------------------------------
// l'(x | mu, sigma) = l(x | mu, sigma) * Gauss(mu_obs | mu, 0.2)
// event observables
RooRealVar x("x", "x", -10, 10);
// parameters
RooRealVar mu("mu", "mu", 0.0, -10, 10);
RooRealVar sigma("sigma", "sigma", 1.0, 0.1, 2.0);
// Gaussian model for event observables
RooGaussian gauss("gauss", "gauss", x, mu, sigma);
// global observables (which are not parameters so they are constant)
RooRealVar mu_obs("mu_obs", "mu_obs", 1.0, -10, 10);
mu_obs.setConstant();
// note: alternatively, one can create a constant with default limits using RooRealVar("mu_obs", "mu_obs", 1.0)
// constraint pdf
RooGaussian constraint("constraint", "constraint", mu_obs, mu, 0.1);
// full pdf including constraint pdf
RooProdPdf model("model", "model", {gauss, constraint});
// Generating toy data with randomized global observables
// ------------------------------------------------------
// For most toy-based statistical procedures, it is necessary to also
// randomize the global observable when generating toy datasets.
//
// To that end, let's generate a single event from the model and take the
// global observable value (the same is done in the RooStats:ToyMCSampler
// class):
std::unique_ptr<RooDataSet> dataGlob{model.generate({mu_obs}, 1)};
// Next, we temporarily set the value of mu_obs to the randomized value for
// generating our toy dataset:
double mu_obs_orig_val = mu_obs.getVal();
RooArgSet{mu_obs}.assign(*dataGlob->get(0));
// Actually generate the toy dataset. We don't generate too many events,
// otherwise, the constraint will not have much weight in the fit and the
// result looks like it's unaffected by it.
std::unique_ptr<RooDataSet> data{model.generate({x}, 50)};
// When fitting the toy dataset, it is important to set the global
// observables in the fit to the values that were used to generate the toy
// dataset. To facilitate the bookkeeping of global observable values, you
// can attach a snapshot with the current global observable values to the
// dataset like this (new feature introduced in ROOT 6.26):
data->setGlobalObservables({mu_obs});
// reset original mu_obs value
mu_obs.setVal(mu_obs_orig_val);
// Fitting a model with global observables
// ---------------------------------------
// Create snapshot of original parameters to reset parameters after fitting
RooArgSet modelParameters;
model.getParameters(data->get(), modelParameters);
RooArgSet origParameters;
modelParameters.snapshot(origParameters);
using FitRes = std::unique_ptr<RooFitResult>;
// When you fit a model that includes global observables, you need to
// specify them in the call to RooAbsPdf::fitTo with the
// RooFit::GlobalObservables command argument. By default, the global
// observable values attached to the dataset will be prioritized over the
// values in the model, so the following fit correctly uses the randomized
// global observable values from the toy dataset:
std::cout << "1. model.fitTo(*data, GlobalObservables(mu_obs))\n";
std::cout << "------------------------------------------------\n";
FitRes res1{model.fitTo(*data, GlobalObservables(mu_obs), PrintLevel(-1), Save())};
res1->Print();
modelParameters.assign(origParameters);
// In our example, the set of global observables is attached to the toy
// dataset. In this case, you can actually drop the GlobalObservables()
// command argument, because the global observables are automatically
// figured out from the data set (this fit result should be identical to the
// previous one).
std::cout << "2. model.fitTo(*data)\n";
std::cout << "---------------------\n";
FitRes res2{model.fitTo(*data, PrintLevel(-1), Save())};
res2->Print();
modelParameters.assign(origParameters);
// If you want to explicitly ignore the global observables in the dataset,
// you can do that by specifying GlobalObservablesSource("model"). Keep in
// mind that now it's also again your responsability to define the set of
// global observables.
std::cout << "3. model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource(\"model\"))\n";
std::cout << "------------------------------------------------\n";
FitRes res3{model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource("model"), PrintLevel(-1), Save())};
res3->Print();
modelParameters.assign(origParameters);
}
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.
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition: RooAbsPdf.h:60
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:179
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
static RooMsgService & instance()
Return reference to singleton instance.
StreamConfig & getStream(Int_t id)
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:33
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
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: Common.h:18
@ Minimization
Definition: RooGlobalFunc.h:62
static constexpr double gauss
void removeTopic(RooFit::MsgTopic oldTopic)
1. model.fitTo(*data, GlobalObservables(mu_obs))
------------------------------------------------
RooFitResult: minimized FCN value: 68.2482, estimated distance to minimum: 1.96065e-06
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu 5.2717e-02 +/- 8.11e-02
sigma 9.7190e-01 +/- 9.73e-02
2. model.fitTo(*data)
---------------------
RooFitResult: minimized FCN value: 68.2482, estimated distance to minimum: 1.96065e-06
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu 5.2717e-02 +/- 8.11e-02
sigma 9.7190e-01 +/- 9.73e-02
3. model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource("model"))
------------------------------------------------
RooFitResult: minimized FCN value: 83.7181, estimated distance to minimum: 1.33582e-06
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu 7.4744e-01 +/- 9.68e-02
sigma 1.2451e+00 +/- 1.38e-01
Date
January 2022

Definition in file rf613_global_observables.C.