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 Tuesday, March 19, 2024 at 07:17 PM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.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;
Observables
RooRealVar t("t", "time", -1., 15.);
RooRealVar cosa("cosa", "cos(alpha)", -1., 1.);
RooRealVar tau("tau", "#tau", 1.5);
RooRealVar deltaGamma("deltaGamma", "deltaGamma", 0.3);
RooFormulaVar coshG("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", {t, tau, deltaGamma});
RooFormulaVar sinhG("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", {t, tau, deltaGamma});
Construct polynomial amplitudes in cos(a)
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)
RooProduct ampl1("ampl1", "amplitude 1", {poly1, coshG});
RooProduct ampl2("ampl2", "amplitude 2", {poly2, sinhG});
Amplitude strengths
RooRealVar f1("f1", "f1", 1, 0, 2);
RooRealVar f2("f2", "f2", 0.5, 0, 2);
Construct pdf
RooRealSumPdf pdf("pdf", "pdf", RooArgList(ampl1, ampl2), RooArgList(f1, f2));
Generate some toy data from pdf
std::unique_ptr<RooDataSet> data{pdf.generate({t, cosa}, 10000)};
[#1] INFO:NumericIntegration -- RooRealIntegral::init(coshGBasis_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t) [#1] INFO:NumericIntegration -- RooRealIntegral::init(sinhGBasis_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t) [#1] INFO:NumericIntegration -- RooRealIntegral::init(coshGBasis_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t) [#1] INFO:NumericIntegration -- RooRealIntegral::init(sinhGBasis_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
input_line_59:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration std::unique_ptr<RooDataSet> data{pdf.generate({t, cosa}, 10000)}; ^
Fit pdf to toy data with only amplitude strength floating
pdf.fitTo(*data, PrintLevel(-1));
input_line_60:2:13: error: reference to 'data' is ambiguous pdf.fitTo(*data, PrintLevel(-1)); ^ input_line_59:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{pdf.generate({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 2D plots of amplitudes
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);
Make projection on t, plot data, pdf and its components Note component projections may be larger than sum because amplitudes can be negative
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_62:3:1: error: reference to 'data' is ambiguous data->plotOn(frame1); ^ input_line_59:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{pdf.generate({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
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_63:3:1: error: reference to 'data' is ambiguous data->plotOn(frame2); ^ input_line_59:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{pdf.generate({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_63:12:1: error: use of undeclared identifier 'frame1' frame1->GetYaxis()->SetTitleOffset(1.4); ^ input_line_63:13:1: error: use of undeclared identifier 'frame1' frame1->Draw(); ^
Draw all canvases
gROOT->GetListOfCanvases()->Draw()