rf705_linearmorph

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, November 30, 2022 at 11:24 AM.

In [1]:
%%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;

Create end point pdf shapes

Observable

In [2]:
RooRealVar x("x", "x", -20, 20);

Lower end point shape: a Gaussian

In [3]:
RooRealVar g1mean("g1mean", "g1mean", -10);
RooGaussian g1("g1", "g1", x, g1mean, RooConst(2));

Upper end point shape: a Polynomial

In [4]:
RooPolynomial g2("g2", "g2", x, RooArgSet(-0.03, -0.001));

Create interpolating pdf

Create interpolation variable

In [5]:
RooRealVar alpha("alpha", "alpha", 0, 1.0);

Specify sampling density on observable and interpolation variable

In [6]:
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

In [7]:
RooIntegralMorph lmorph("lmorph", "lmorph", g1, g2, x, alpha);

Plot interpolating pdf at various alpha

Show end points as blue curves

In [8]:
RooPlot *frame1 = x.frame();
g1.plotOn(frame1);
g2.plotOn(frame1);

Show interpolated shapes in red

In [9]:
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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 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 0x7f9aa52fc9d0 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0

Show 2D distribution of pdf(x,alpha)

Create 2D histogram

In [10]:
TH1 *hh = lmorph.createHistogram("hh", x, Binning(40), YVar(alpha, Binning(40)));
hh->SetLineColor(kBlue);
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f9aa5328c20 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x_alpha for nset (x,alpha) with code 0 from preexisting content.

Fit pdf to dataset with alpha=0.8

Generate a toy dataset at alpha = 0.8

In [11]:
alpha = 0.8;
RooDataSet *data = lmorph.generate(x, 1000);
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(lmorph) creating new cache 0x7f9aa54e41a0 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 0x7f9aa549d270 with pdf g1_MORPH_g2_CACHE_Obs[x]_NORM_x for nset (x) with code 0 from preexisting content.
input_line_59:3:1: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
RooDataSet *data = lmorph.generate(x, 1000);
^

Fit pdf to toy data

In [12]:
lmorph.setCacheAlpha(true);
lmorph.fitTo(*data, Verbose(true));
input_line_60:3:15: error: reference to 'data' is ambiguous
lmorph.fitTo(*data, Verbose(true));
              ^
input_line_59:3:13: note: candidate found by name lookup is '__cling_N530::data'
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

In [13]:
RooPlot *frame2 = x.frame(Bins(100));
data->plotOn(frame2);
lmorph.plotOn(frame2);
input_line_61:3:1: error: reference to 'data' is ambiguous
data->plotOn(frame2);
^
input_line_59:3:13: note: candidate found by name lookup is '__cling_N530::data'
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
    ^

Scan -log(L) vs alpha

Show scan -log(L) of dataset w.r.t alpha

In [14]:
RooPlot *frame3 = alpha.frame(Bins(100), Range(0.1, 0.9));

Make 2D pdf of histogram

In [15]:
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_63:2:52: error: reference to 'data' is ambiguous
 std::unique_ptr<RooAbsReal> nll{lmorph.createNLL(*data)};
                                                   ^
input_line_59:3:13: note: candidate found by name lookup is '__cling_N530::data'
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_63:23:1: error: use of undeclared identifier 'frame2'
frame2->GetYaxis()->SetTitleOffset(1.4);
^
input_line_63:24:1: error: use of undeclared identifier 'frame2'
frame2->Draw();
^

Draw all canvases

In [16]:
gROOT->GetListOfCanvases()->Draw()