Setting up an extended maximum likelihood fit.
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:17 AM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit;
Declare observable x
RooRealVar x("x", "x", 0, 10);
Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean", "mean of gaussians", 5);
RooRealVar sigma1("sigma1", "width of gaussians", 0.5);
RooRealVar sigma2("sigma2", "width of gaussians", 1);
RooGaussian sig1("sig1", "Signal component 1", x, mean, sigma1);
RooGaussian sig2("sig2", "Signal component 2", x, mean, sigma2);
[#0] WARNING:InputArguments -- The parameter 'sigma1' with range [-inf, inf] of the RooGaussian 'sig1' exceeds the safe range of (0, inf). Advise to limit its range. [#0] WARNING:InputArguments -- The parameter 'sigma2' with range [-inf, inf] of the RooGaussian 'sig2' exceeds the safe range of (0, inf). Advise to limit its range.
Build Chebychev polynomial pdf
RooRealVar a0("a0", "a0", 0.5, 0., 1.);
RooRealVar a1("a1", "a1", 0.2, 0., 1.);
RooChebychev bkg("bkg", "Background", x, RooArgSet(a0, a1));
Sum the signal components into a composite signal pdf
RooRealVar sig1frac("sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.);
RooAddPdf sig("sig", "Signal", RooArgList(sig1, sig2), sig1frac);
Sum the composite signal and background into an extended pdf nsigsig+nbkgbkg
RooRealVar nsig("nsig", "number of signal events", 500, 0., 10000);
RooRealVar nbkg("nbkg", "number of background events", 500, 0, 10000);
RooAddPdf model("model", "(g1+g2)+a", RooArgList(bkg, sig), RooArgList(nbkg, nsig));
Generate a data sample of expected number events in x from model = model.expectedEvents() = nsig+nbkg
std::unique_ptr<RooDataSet> data{model.generate(x)};
input_line_53: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{model.generate(x)}; ^
Fit model to data, extended ML term automatically included
model.fitTo(*data, PrintLevel(-1));
input_line_54:2:15: error: reference to 'data' is ambiguous model.fitTo(*data, PrintLevel(-1)); ^ input_line_53:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.generate(x)}; ^ /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 data and PDF overlaid, use expected number of events for pdf projection normalization rather than observed number of events (==data->numEntries())
RooPlot *xframe = x.frame(Title("extended ML fit example"));
data->plotOn(xframe);
model.plotOn(xframe, Normalization(1.0, RooAbsReal::RelativeExpected));
input_line_55:3:1: error: reference to 'data' is ambiguous data->plotOn(xframe); ^ input_line_53:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.generate(x)}; ^ /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 ^
Overlay the background component of model with a dashed line
model.plotOn(xframe, Components(bkg), LineStyle(kDashed), Normalization(1.0, RooAbsReal::RelativeExpected));
input_line_56:2:50: error: cannot take the address of an rvalue of type 'ELineStyle' model.plotOn(xframe, Components(bkg), LineStyle(kDashed), Normalization(1.0, RooAbsReal::RelativeExpected)); ^~~~~~~ Error while creating dynamic expression for: model.plotOn(xframe, Components(bkg), LineStyle(kDashed), Normalization(1., RooAbsReal::RelativeExpected))
Overlay the background+sig2 components of model with a dotted line
model.plotOn(xframe, Components(RooArgSet(bkg, sig2)), LineStyle(kDotted),
Normalization(1.0, RooAbsReal::RelativeExpected));
input_line_57:2:67: error: cannot take the address of an rvalue of type 'ELineStyle' model.plotOn(xframe, Components(RooArgSet(bkg, sig2)), LineStyle(kDotted), ^~~~~~~ Error while creating dynamic expression for: model.plotOn(xframe, Components(RooArgSet(bkg, sig2)), LineStyle(kDotted), Normalization(1., RooAbsReal::RelativeExpected))
Print structure of composite pdf
model.Print("t");
0x7fb7323727d0 RooAddPdf::model = 0.9/1 [Auto,Clean] 0x7fb732a4b7d0/V- RooChebychev::bkg = 0.8 [Auto,Dirty] 0x7fb748009000/V- RooRealVar::x = 5 0x7fb732a4b000/V- RooRealVar::a0 = 0.5 0x7fb732a4b3e8/V- RooRealVar::a1 = 0.2 0x7fb7323723e8/V- RooRealVar::nbkg = 500 0x7fb7323753e8/V- RooAddPdf::sig = 1/1 [Auto,Clean] 0x7fb732a4ebb8/V- RooGaussian::sig1 = 1 [Auto,Dirty] 0x7fb748009000/V- RooRealVar::x = 5 0x7fb732a4e000/V- RooRealVar::mean = 5 0x7fb732a4e3e8/V- RooRealVar::sigma1 = 0.5 0x7fb732375000/V- RooRealVar::sig1frac = 0.8 0x7fb732a4f110/V- RooGaussian::sig2 = 1 [Auto,Dirty] 0x7fb748009000/V- RooRealVar::x = 5 0x7fb732a4e000/V- RooRealVar::mean = 5 0x7fb732a4e7d0/V- RooRealVar::sigma2 = 1 0x7fb732372000/V- RooRealVar::nsig = 500
Associated nsig/nbkg as expected number of events with sig/bkg
RooExtendPdf esig("esig", "extended signal pdf", sig, nsig);
RooExtendPdf ebkg("ebkg", "extended background pdf", bkg, nbkg);
Construct sum of two extended pdf (no coefficients required)
RooAddPdf model2("model2", "(g1+g2)+a", RooArgList(ebkg, esig));
Draw the frame on the canvas
new TCanvas("rf202_composite", "rf202_composite", 600, 600);
gPad->SetLeftMargin(0.15);
xframe->GetYaxis()->SetTitleOffset(1.4);
xframe->Draw();
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
%jsroot on
gROOT->GetListOfCanvases()->Draw()