Special pdf's: linear interpolation between pdf shapes using the 'Alex Read' algorithm
Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, April 17, 2024 at 11:19 AM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooIntegralMorph.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "TH1.h"
using namespace RooFit;
Observable
RooRealVar x("x", "x", -20, 20);
Lower end point shape: a Gaussian
RooRealVar g1mean("g1mean", "g1mean", -10);
RooGaussian g1("g1", "g1", x, g1mean, 2.0);
Upper end point shape: a Polynomial
RooPolynomial g2("g2", "g2", x, RooArgSet(-0.03, -0.001));
Create interpolation variable
RooRealVar alpha("alpha", "alpha", 0, 1.0);
Specify sampling density on observable and interpolation variable
x.setBins(1000, "cache");
alpha.setBins(50, "cache");
Construct interpolating pdf in (x,a) represent g1(x) at a=a_min and g2(x) at a=a_max
RooIntegralMorph lmorph("lmorph", "lmorph", g1, g2, x, alpha);
Show end points as blue curves
RooPlot *frame1 = x.frame();
g1.plotOn(frame1);
g2.plotOn(frame1);
Show interpolated shapes in red
alpha.setVal(0.125);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.25);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.375);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.50);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.625);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.75);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.875);
lmorph.plotOn(frame1, LineColor(kRed));
alpha.setVal(0.95);
lmorph.plotOn(frame1, LineColor(kRed));
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99a8a20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0
Create 2D histogram
TH1 *hh = lmorph.createHistogram("hh", x, Binning(40), YVar(alpha, Binning(40)));
hh->SetLineColor(kBlue);
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f9a68810 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x_alpha for nset (x,alpha) with code 0 from preexisting content.
Generate a toy dataset at alpha = 0.8
alpha = 0.8;
std::unique_ptr<RooDataSet> data{lmorph.generate(x, 1000)};
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f9a991d0 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f20f99cc6d0 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 from preexisting content.
input_line_57:3:1: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration std::unique_ptr<RooDataSet> data{lmorph.generate(x, 1000)}; ^
Fit pdf to toy data
lmorph.setCacheAlpha(true);
lmorph.fitTo(*data, Verbose(true), PrintLevel(-1));
input_line_58:3:15: error: reference to 'data' is ambiguous lmorph.fitTo(*data, Verbose(true), PrintLevel(-1)); ^ input_line_57:3:29: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{lmorph.generate(x, 1000)}; ^ /usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data' data(initializer_list<_Tp> __il) noexcept ^ /usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data' data(_Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data' data(const _Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data' data(_Tp (&__array)[_Nm]) noexcept ^
Plot fitted pdf and data overlaid
RooPlot *frame2 = x.frame(Bins(100));
data->plotOn(frame2);
lmorph.plotOn(frame2);
input_line_59:3:1: error: reference to 'data' is ambiguous data->plotOn(frame2); ^ input_line_57:3:29: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{lmorph.generate(x, 1000)}; ^ /usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data' data(initializer_list<_Tp> __il) noexcept ^ /usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data' data(_Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data' data(const _Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data' data(_Tp (&__array)[_Nm]) noexcept ^
Show scan -log(L) of dataset w.r.t alpha
RooPlot *frame3 = alpha.frame(Bins(100), Range(0.1, 0.9));
Make 2D pdf of histogram
std::unique_ptr<RooAbsReal> nll{lmorph.createNLL(*data)};
nll->plotOn(frame3, ShiftToZero());
lmorph.setCacheAlpha(false);
TCanvas *c = new TCanvas("rf705_linearmorph", "rf705_linearmorph", 800, 800);
c->Divide(2, 2);
c->cd(1);
gPad->SetLeftMargin(0.15);
frame1->GetYaxis()->SetTitleOffset(1.6);
frame1->Draw();
c->cd(2);
gPad->SetLeftMargin(0.20);
hh->GetZaxis()->SetTitleOffset(2.5);
hh->Draw("surf");
c->cd(3);
gPad->SetLeftMargin(0.15);
frame3->GetYaxis()->SetTitleOffset(1.4);
frame3->Draw();
c->cd(4);
gPad->SetLeftMargin(0.15);
frame2->GetYaxis()->SetTitleOffset(1.4);
frame2->Draw();
return;
input_line_61:2:52: error: reference to 'data' is ambiguous std::unique_ptr<RooAbsReal> nll{lmorph.createNLL(*data)}; ^ input_line_57:3:29: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{lmorph.generate(x, 1000)}; ^ /usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data' data(initializer_list<_Tp> __il) noexcept ^ /usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data' data(_Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data' data(const _Container& __cont) noexcept(noexcept(__cont.data())) ^ /usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data' data(_Tp (&__array)[_Nm]) noexcept ^ input_line_61:23:1: error: use of undeclared identifier 'frame2' frame2->GetYaxis()->SetTitleOffset(1.4); ^ input_line_61:24:1: error: use of undeclared identifier 'frame2' frame2->Draw(); ^
Draw all canvases
gROOT->GetListOfCanvases()->Draw()