Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf616_morphing.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4/// Use Morphing in RooFit.
5///
6/// This tutorial shows how to use template morphing inside RooFit. As input we have several
7/// Gaussian distributions. The output is one gaussian, with a specific mean value.
8/// Since likelihoods are often used within the framework of morphing, we provide a
9/// way to estimate the negative log likelihood (nll).
10///
11/// Based on example of Kyle Cranmer https://gist.github.com/cranmer/b67830e46d53d5f7cf2d.
12///
13/// \macro_image
14/// \macro_code
15/// \macro_output
16///
17/// \date August 2024
18/// \author Robin Syring
19
20#include "RooRealVar.h"
21#include "RooRealVar.h"
22#include "RooWorkspace.h"
23#include "RooGaussian.h"
24#include "RooUniform.h"
25#include "RooDataSet.h"
26#include "RooPlot.h"
28#include "RooAbsPdf.h"
29
30using namespace RooFit;
31
32// Number of samples to fill the histograms
33const int n_samples = 1000;
34
35// Kills warning massages
37
38// Define the morphing routine
39RooPlot *perform_morphing(RooWorkspace &ws, RooMomentMorphFuncND::Setting setting, double sigma)
40{
41 // Get Variables from the workspace
42 RooRealVar *x_var = ws.var("x");
43 RooRealVar *mu_var = ws.var("mu");
44 RooAbsPdf *gauss = ws.pdf("gauss");
45
46 // Initialize a plot
47 RooPlot *frame1 = x_var->frame();
48
49 // Define binning for morphing
50 RooMomentMorphFuncND::Grid grid(RooBinning(4, 0.0, 4.0));
51
52 // Set binning of histograms, has to be customized for optimal results
53 x_var->setBins(50);
54
55 std::vector<int> parampoints = {0, 1, 2, 3, 4};
56
57 for (auto i : parampoints) {
58 // Define the sampled gaussians
59 RooRealVar mu_help(Form("mu%d", i), Form("mu%d", i), i);
60 // Use * because RooGaussian expects objects no pointers
61 RooGaussian help(Form("g%d", i), Form("g%d", i), *x_var, mu_help, sigma);
62 ws.import(help, Silence(true));
63
64 // Fill the histograms use a unique pointer to prevent memory leaks
65 std::unique_ptr<RooDataHist> hist1{
66 dynamic_cast<RooDataHist *>(ws.pdf(Form("g%d", i))->generateBinned(*x_var, 100 * n_samples))};
67
68 // Add the value 1 to each bin
69 for (int i_bin = 0; i_bin < hist1->numEntries(); ++i_bin) {
70 const RooArgSet *binContent = hist1->get(i_bin);
71 hist1->add(*binContent, 1.0);
72 }
73
74 // Add the pdf to the workspace, the inOrder of 1 is necessary for calculation of the nll
75 // Adjust it to 0 to see binning
76 ws.import(RooHistPdf(Form("histpdf%d", i), Form("histpdf%d", i), *x_var, *hist1, 1), Silence(true));
77
78 // Plot and add the pdf to the grid
79 RooAbsPdf *pdf = ws.pdf(Form("histpdf%d", i));
80 pdf->plotOn(frame1);
81 grid.addPdf(*pdf, i);
82 }
83
84 // Create the morphing
85 RooMomentMorphFuncND morph_func("morpf_func", "morph_func", RooArgList(*mu_var), RooArgList(*x_var), grid, setting);
86
87 // Normalizing the morphed object to be a pdf, set it false to prevent warning messages and gain computational speed
88 // up
89 morph_func.setPdfMode();
90
91 // Creating the morphed pdf
92 RooWrapperPdf morph("morph", "morph", morph_func, true);
93 ws.import(morph, Silence(true));
94 RooAbsPdf *morph_ = ws.pdf("morph");
95 morph_->plotOn(frame1, LineColor(kRed));
96
97 return frame1;
98}
99
100// Define the workspace
101std::unique_ptr<RooWorkspace> build_ws(double mu_observed, double sigma)
102{
103 auto ws = std::make_unique<RooWorkspace>();
104 ws->factory(Form("Gaussian::gauss(x[-5,15],mu[%f,0,4], %f)", mu_observed, sigma));
105 return ws;
106}
107
108// Do the example
109void rf616_morphing()
110{
111 // Define the 'observed' mu
112 double mu_observed = 2.5;
113 double sigma = 1.5;
114
115 // Import variables from workspace
116 std::unique_ptr<RooWorkspace> ws = build_ws(mu_observed, sigma);
117
118 RooPlot *frame1 = perform_morphing(*ws, RooMomentMorphFuncND::Linear, sigma);
119
120 RooRealVar *x_var = ws->var("x");
121 RooRealVar *mu_var = ws->var("mu");
122 RooAbsPdf *gauss = ws->pdf("gauss");
123 RooDataSet *obs_data = gauss->generate(*x_var, n_samples);
124
125 // Create the exact negative log likelihood function for Gaussian model
126 RooAbsReal *nll_gauss = gauss->createNLL(*obs_data);
127
128 // Create the morphed negative log likelihood function
129 // TODO: Fix RooAddPdf::fixCoefNormalization(nset) warnings with new CPU backend
130 RooAbsReal *nll_morph = ws->pdf("morph")->createNLL(*obs_data, RooFit::EvalBackend("legacy"));
131
132 // Plot the negative logarithmic summed likelihood
133 RooPlot *frame2 = mu_var->frame(Title("Negative log Likelihood"));
134 nll_gauss->plotOn(frame2, LineColor(kBlue), ShiftToZero(), Name("gauss"));
135 nll_morph->plotOn(frame2, LineColor(kRed), ShiftToZero(), Name("morph"));
136
137 TCanvas *c = new TCanvas("rf616_morphing", "rf616_morphing", 800, 400);
138 c->Divide(2);
139 c->cd(1);
140 gPad->SetLeftMargin(0.15);
141 frame1->GetYaxis()->SetTitleOffset(1.4);
142 frame1->Draw();
143 c->cd(2);
144 gPad->SetLeftMargin(0.15);
145 frame2->GetYaxis()->SetTitleOffset(1.8);
146 frame2->Draw();
147
148 // Compute the minimum of the nll via minuit
149 std::vector<RooAbsReal *> nlls = {nll_gauss, nll_morph};
150 for (auto nll : nlls) {
151 RooMinimizer minimizer(*nll);
152 minimizer.setPrintLevel(-1);
153 minimizer.minimize("Minuit2");
154 RooFitResult *result = minimizer.save();
155 result->Print();
156 }
157}
#define c(i)
Definition RSha256.hxx:101
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
#define gPad
Storage_t const & get() const
Const access to the underlying stl container.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:163
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:124
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
virtual RooFit::OwningPtr< RooDataHist > generateBinned(const RooArgSet &whatVars, double nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}) const
As RooAbsPdf::generateBinned(const RooArgSet&, const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,...
Definition RooAbsPdf.h:110
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const
Plot (project) PDF on specified frame.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
Implements a RooAbsBinning in terms of an array of boundary values, posing no constraints on the choi...
Definition RooBinning.h:27
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
Container class to hold unbinned data.
Definition RooDataSet.h:34
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
A propability density function sampled from a multidimensional histogram.
Definition RooHistPdf.h:30
Wrapper class around ROOT::Math::Minimizer that provides a seamless interface between the minimizer f...
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
TAxis * GetYaxis() const
Definition RooPlot.cxx:1224
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:597
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
Persistable container for RooFit projects.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
RooFactoryWSTool & factory()
Return instance to factory tool.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
The RooWrapperPdf is a class that can be used to convert a function into a PDF.
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
The Canvas class.
Definition TCanvas.h:23
RooCmdArg Silence(bool flag=true)
RooCmdArg ShiftToZero()
RooCmdArg LineColor(Color_t color)
const Double_t sigma
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:64
const char * Name
Definition TXMLSetup.cxx:67
const char * Title
Definition TXMLSetup.cxx:68