rf704_amplitudefit

Special pdf's: using a pdf defined by a sum of real-valued amplitude components

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 "RooTruthModel.h"
#include "RooFormulaVar.h"
#include "RooRealSumPdf.h"
#include "RooPolyVar.h"
#include "RooProduct.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit;

Setup 2D amplitude functions

Observables

In [2]:
RooRealVar t("t", "time", -1., 15.);
RooRealVar cosa("cosa", "cos(alpha)", -1., 1.);

Use RooTruthModel to obtain compiled implementation of sinh/cosh modulated decay functions

In [3]:
RooRealVar tau("tau", "#tau", 1.5);
RooRealVar deltaGamma("deltaGamma", "deltaGamma", 0.3);
RooTruthModel truthModel("tm", "tm", t);
RooFormulaVar coshGBasis("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", RooArgList(t, tau, deltaGamma));
RooFormulaVar sinhGBasis("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", RooArgList(t, tau, deltaGamma));
RooAbsReal *coshGConv = truthModel.convolution(&coshGBasis, &t);
RooAbsReal *sinhGConv = truthModel.convolution(&sinhGBasis, &t);

Construct polynomial amplitudes in cos(a)

In [4]:
RooPolyVar poly1("poly1", "poly1", cosa, RooArgList(0.5, 0.2, 0.2), 0);
RooPolyVar poly2("poly2", "poly2", cosa, RooArgList(1.0, -0.2, 3.0), 0);

Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)

In [5]:
RooProduct ampl1("ampl1", "amplitude 1", RooArgSet(poly1, *coshGConv));
RooProduct ampl2("ampl2", "amplitude 2", RooArgSet(poly2, *sinhGConv));

Construct amplitude sum pdf

Amplitude strengths

In [6]:
RooRealVar f1("f1", "f1", 1, 0, 2);
RooRealVar f2("f2", "f2", 0.5, 0, 2);

Construct pdf

In [7]:
RooRealSumPdf pdf("pdf", "pdf", RooArgList(ampl1, ampl2), RooArgList(f1, f2));

Generate some toy data from pdf

In [8]:
RooDataSet *data = pdf.generate(RooArgSet(t, cosa), 10000);
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
input_line_62:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooDataSet *data = pdf.generate(RooArgSet(t, cosa), 10000);
 ^

Fit pdf to toy data with only amplitude strength floating

In [9]:
pdf.fitTo(*data);
input_line_63:2:13: error: reference to 'data' is ambiguous
 pdf.fitTo(*data);
            ^
input_line_62:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = pdf.generate(RooArgSet(t, cosa), 10000);
             ^
/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 amplitude sum pdf

Make 2D plots of amplitudes

In [10]:
TH1 *hh_cos = ampl1.createHistogram("hh_cos", t, Binning(50), YVar(cosa, Binning(50)));
TH1 *hh_sin = ampl2.createHistogram("hh_sin", t, Binning(50), YVar(cosa, Binning(50)));
hh_cos->SetLineColor(kBlue);
hh_sin->SetLineColor(kRed);
[#0] WARNING:InputArguments -- RooAbsReal::createHistogram(ampl1) WARNING extended mode requested for a non-pdf object, ignored
[#0] WARNING:InputArguments -- RooAbsReal::createHistogram(ampl2) WARNING extended mode requested for a non-pdf object, ignored

Make projection on t, plot data, pdf and its components Note component projections may be larger than sum because amplitudes can be negative

In [11]:
RooPlot *frame1 = t.frame();
data->plotOn(frame1);
pdf.plotOn(frame1);
pdf.plotOn(frame1, Components(ampl1), LineStyle(kDashed));
pdf.plotOn(frame1, Components(ampl2), LineStyle(kDashed), LineColor(kRed));
input_line_65:3:1: error: reference to 'data' is ambiguous
data->plotOn(frame1);
^
input_line_62:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = pdf.generate(RooArgSet(t, cosa), 10000);
             ^
/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
    ^

Make projection on cosa, plot data, pdf and its components Note that components projection may be larger than sum because amplitudes can be negative

In [12]:
RooPlot *frame2 = cosa.frame();
data->plotOn(frame2);
pdf.plotOn(frame2);
pdf.plotOn(frame2, Components(ampl1), LineStyle(kDashed));
pdf.plotOn(frame2, Components(ampl2), LineStyle(kDashed), LineColor(kRed));

TCanvas *c = new TCanvas("rf704_amplitudefit", "rf704_amplitudefit", 800, 800);
c->Divide(2, 2);
c->cd(1);
gPad->SetLeftMargin(0.15);
frame1->GetYaxis()->SetTitleOffset(1.4);
frame1->Draw();
c->cd(2);
gPad->SetLeftMargin(0.15);
frame2->GetYaxis()->SetTitleOffset(1.4);
frame2->Draw();
c->cd(3);
gPad->SetLeftMargin(0.20);
hh_cos->GetZaxis()->SetTitleOffset(2.3);
hh_cos->Draw("surf");
c->cd(4);
gPad->SetLeftMargin(0.20);
hh_sin->GetZaxis()->SetTitleOffset(2.3);
hh_sin->Draw("surf");
input_line_66:3:1: error: reference to 'data' is ambiguous
data->plotOn(frame2);
^
input_line_62:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = pdf.generate(RooArgSet(t, cosa), 10000);
             ^
/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_66:12:1: error: use of undeclared identifier 'frame1'
frame1->GetYaxis()->SetTitleOffset(1.4);
^
input_line_66:13:1: error: use of undeclared identifier 'frame1'
frame1->Draw();
^

Draw all canvases

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