Extended maximum likelihood fit in multiple ranges.
When an extended pdf and multiple ranges are used, the RooExtendPdf cannot correctly interpret the coefficients used for extension. This can be solved by using a RooAddPdf for extending the model.
Author: Stephan Hageboeck
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:15 PM.
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit ;
Declare observable x
RooRealVar x("x","x",0,11) ;
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) ;
Define signal range in which events counts are to be defined
x.setRange("signalRange",4,6) ;
[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'signalRange' created with bounds [4,6]
Associated nsig/nbkg as expected number of events with sig/bkg in_the_range "signalRange"
RooRealVar nsig("nsig","number of signal events in signalRange",500,0.,10000) ;
RooRealVar nbkg("nbkg","number of background events in signalRange",500,0,10000) ;
Use AddPdf to extend the model. Giving as many coefficients as pdfs switches on extension.
RooAddPdf model("model","(g1+g2)+a", RooArgList(bkg,sig), RooArgList(nbkg,nsig)) ;
Generate 1000 events from model so that nsig,nbkg come out to numbers <<500 in fit
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
auto canv = new TCanvas("Canvas", "Canvas", 1500, 600);
canv->Divide(3,1);
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{model.generate(x,1000)}; ^
canv->cd(1);
Perform unbinned ML fit to data, full range
IMPORTANT: The model needs to be copied when fitting with different ranges because the interpretation of the coefficients is tied to the fit range that's used in the first fit
RooAddPdf model1(model);
std::unique_ptr<RooFitResult> r{model1.fitTo(*data,Save(), PrintLevel(-1))};
r->Print() ;
RooPlot * frame = x.frame(Title("Full range fitted"));
data->plotOn(frame);
model1.plotOn(frame, VisualizeError(*r));
model1.plotOn(frame);
model1.paramOn(frame);
frame->Draw();
input_line_59:3:47: error: reference to 'data' is ambiguous std::unique_ptr<RooFitResult> r{model1.fitTo(*data,Save(), PrintLevel(-1))}; ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.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_59:7:1: error: reference to 'data' is ambiguous data->plotOn(frame); ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.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 ^
canv->cd(2);
x.setRange("left", 0., 4.);
x.setRange("right", 6., 10.);
RooAddPdf model2(model);
std::unique_ptr<RooFitResult> r2{model2.fitTo(*data, Range("left,right"), Save(), PrintLevel(-1))};
r2->Print();
RooPlot * frame2 = x.frame(Title("Fit in left/right sideband"));
data->plotOn(frame2);
model2.plotOn(frame2, VisualizeError(*r2));
model2.plotOn(frame2);
model2.paramOn(frame2);
frame2->Draw();
input_line_60:7:48: error: reference to 'data' is ambiguous std::unique_ptr<RooFitResult> r2{model2.fitTo(*data, Range("left,right"), Save(), PrintLevel(-1))}; ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.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_60:12: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{model.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 ^
Note how restricting the region to only the left tail increases the fit uncertainty
canv->cd(3);
x.setRange("leftToMiddle", 0., 5.);
RooAddPdf model3(model);
std::unique_ptr<RooFitResult> r3{model3.fitTo(*data, Range("leftToMiddle"), Save(), PrintLevel(-1))};
r3->Print();
RooPlot * frame3 = x.frame(Title("Fit from left to middle"));
data->plotOn(frame3);
model3.plotOn(frame3, VisualizeError(*r3));
model3.plotOn(frame3);
model3.paramOn(frame3);
frame3->Draw();
canv->Draw();
input_line_61:6:48: error: reference to 'data' is ambiguous std::unique_ptr<RooFitResult> r3{model3.fitTo(*data, Range("leftToMiddle"), Save(), PrintLevel(-1))}; ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.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_61:11:1: error: reference to 'data' is ambiguous data->plotOn(frame3); ^ input_line_55:2:30: note: candidate found by name lookup is 'data' std::unique_ptr<RooDataSet> data{model.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 ^
Draw all canvases
%jsroot on
gROOT->GetListOfCanvases()->Draw()