Multidimensional models: conditional pdf with per-event errors
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:18 AM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooLandau.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
using namespace RooFit;
Observables
RooRealVar dt("dt", "dt", -10, 10);
RooRealVar dterr("dterr", "per-event error on dt", 0.01, 10);
Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
RooRealVar bias("bias", "bias", 0, -10, 10);
RooRealVar sigma("sigma", "per-event error scale factor", 1, 0.1, 10);
RooGaussModel gm("gm1", "gauss model scaled bt per-event error", dt, bias, sigma, dterr);
Construct decay(dt) (x) gauss1(dt|dterr)
RooRealVar tau("tau", "tau", 1.548);
RooDecay decay_gm("decay_gm", "decay", dt, tau, gm, RooDecay::DoubleSided);
Use landau pdf to get somewhat realistic distribution with long tail
RooLandau pdfDtErr("pdfDtErr", "pdfDtErr", dterr, 1.0, 0.25);
std::unique_ptr<RooDataSet> expDataDterr{pdfDtErr.generate(dterr, 10000)};
Specify external dataset with dterr values to use decay_dm as conditional pdf
std::unique_ptr<RooDataSet> data{decay_gm.generate(dt, ProtoData(*expDataDterr))};
input_line_55: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{decay_gm.generate(dt, ProtoData(*expDataDterr))}; ^
Specify dterr as conditional observable
decay_gm.fitTo(*data, ConditionalObservables(dterr), PrintLevel(-1));
input_line_59:2:18: error: reference to 'data' is ambiguous decay_gm.fitTo(*data, ConditionalObservables(dterr), PrintLevel(-1)); ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{decay_gm.generate(dt, ProtoData(*expDataDterr))}; ^ /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 two-dimensional plot of conditional pdf in (dt,dterr)
TH1 *hh_decay = decay_gm.createHistogram("hh_decay", dt, Binning(50), YVar(dterr, Binning(50)));
hh_decay->SetLineColor(kBlue);
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gm1_conv_exp(-abs(@0)/@1)_dt_tau_[decay_gm]_Int[dt,dterr]) using numeric integrator RooIntegrator1D to calculate Int(dterr)
Plot decay_gm(dt|dterr) at various values of dterr
RooPlot *frame = dt.frame(Title("Slices of decay(dt|dterr) at various dterr"));
for (Int_t ibin = 0; ibin < 100; ibin += 20) {
dterr.setBin(ibin);
decay_gm.plotOn(frame, Normalization(5.));
}
Make projection of data an dt
RooPlot *frame2 = dt.frame(Title("Projection of decay(dt|dterr) on dt"));
data->plotOn(frame2);
input_line_62:3:1: error: reference to 'data' is ambiguous data->plotOn(frame2); ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{decay_gm.generate(dt, ProtoData(*expDataDterr))}; ^ /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 of decay(dt|dterr) on dt.
Instead of integrating out dterr, make a weighted average of curves at values dterr_i as given in the external dataset. (The true argument bins the data before projection to speed up the process)
decay_gm.plotOn(frame2, ProjWData(*expDataDterr, true));
input_line_63:2:26: error: cannot initialize an array element of type 'void *' with an rvalue of type 'RooCmdArg (*)(const RooAbsData &, bool)' decay_gm.plotOn(frame2, ProjWData(*expDataDterr, true)); ^~~~~~~~~
Draw all frames on canvas
TCanvas *c = new TCanvas("rf306_condpereventerrors", "rf306_condperventerrors", 1200, 400);
c->Divide(3);
c->cd(1);
gPad->SetLeftMargin(0.20);
hh_decay->GetZaxis()->SetTitleOffset(2.5);
hh_decay->Draw("surf");
c->cd(2);
gPad->SetLeftMargin(0.15);
frame->GetYaxis()->SetTitleOffset(1.6);
frame->Draw();
c->cd(3);
gPad->SetLeftMargin(0.15);
frame2->GetYaxis()->SetTitleOffset(1.6);
frame2->Draw();
[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { _ZN12__cling_N5301cE, _ZN5cling7runtime8internal15DynamicExprInfoC1EPKcPPvb, __cxx_global_var_initcling_module_324_, $.cling-module-324.__inits.0, _ZN12__cling_N53016__cling_Un1Qu329EPv, __orc_init_func.cling-module-324, _ZN5cling7runtime8internal15DynamicExprInfoC2EPKcPPvb, _ZN3TH18GetZaxisEv, _GLOBAL__sub_I_cling_module_324, _ZN7TObjectnwEm }) } IncrementalExecutor::executeFunction: symbol '_ZN5cling7runtime8internal9EvaluateTIvEET_PNS1_15DynamicExprInfoEPN5clang11DeclContextE' unresolved while linking [cling interface function]! You are probably missing the definition of void cling::runtime::internal::EvaluateT<void>(cling::runtime::internal::DynamicExprInfo*, clang::DeclContext*) Maybe you need to load the corresponding shared library?
Draw all canvases
gROOT->GetListOfCanvases()->Draw()
[runStaticInitializersOnce]: Failed to materialize symbols: { (main, { __orc_init_func.cling-module-324 }) }