ROOT   Reference Guide
Searching...
No Matches
rf204b_extendedLikelihood_rangedFit.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -nodraw
4/// This macro demonstrates how to set up a fit in two ranges for plain
5/// likelihoods and extended likelihoods.
6///
7/// ### 1. Shape fits (plain likelihood)
8///
9/// If you fit a non-extended pdf in two ranges, e.g. pdf->fitTo(data,Range("Range1,Range2")),
10/// it will fit the shapes in the two selected ranges and also take into account the relative
11/// predicted yields in those ranges.
12///
13/// This is useful for example to represent a full-range fit, but with a
14/// blinded signal region inside it.
15///
16///
17/// ### 2. Shape+rate fits (extended likelihood)
18///
19/// If your pdf is extended, i.e. measuring both the distribution in the observable as well
20/// as the event count in the fitted region, some intervention is needed to make fits in ranges
21/// work in a way that corresponds to intuition.
22///
23/// If an extended fit is performed in a sub-range, the observed yield is only that of the subrange, hence
24/// the expected event count will converge to a number that is smaller than what's visible in a plot.
25/// In such cases, it is often preferred to interpret the extended term with respect to the full range
26/// that's plotted, i.e., apply a correction to the extended likelihood term in such a way
27/// that the interpretation of the expected event count remains that of the full range. This can
28/// be done by applying a correcion factor (equal to the fraction of the pdf that is contained in the
29/// fitted range) in the Poisson term that represents the extended likelihood term.
30///
31/// If an extended likelihood fit is performed over *two* sub-ranges, this correction is
32/// even more important: without it, each component likelihood would have a different interpretation
33/// of the expected event count (each corresponding to the count in its own region), and a joint
34/// fit of these regions with different interpretations of the same model parameter results
35/// in a number that is not easily interpreted.
36///
37/// If both regions correct their interpretations such that N_expected refers to the full range,
38/// it is interpreted easily, and consistent in both regions.
39///
40/// This requires that the likelihood model is extended using RooAddPdf in the
41/// form SumPdf = Nsig * sigPdf + Nbkg * bkgPdf.
42///
43/// \macro_image
44/// \macro_code
45/// \macro_output
46///
47/// \authors Stephan Hageboeck, Wouter Verkerke
48
49#include "RooRealVar.h"
50#include "RooExponential.h"
51#include "RooGaussian.h"
53#include "RooDataSet.h"
54#include "RooPlot.h"
55#include "RooExtendPdf.h"
56#include "RooFitResult.h"
57
58#include "TCanvas.h"
59
61{
62 using namespace RooFit;
63
64 // PART 1: Background-only fits
65 // ----------------------------
66
67 // Build plain exponential model
68 RooRealVar x("x", "x", 10, 100);
69 RooRealVar alpha("alpha", "alpha", -0.04, -0.1, -0.0);
70 RooExponential model("model", "Exponential model", x, alpha);
71
72 // Define side band regions and full range
73 x.setRange("LEFT",10,20);
74 x.setRange("RIGHT",60,100);
75
76 x.setRange("FULL",10,100);
77
78 std::unique_ptr<RooDataSet> data{model.generate(x, 10000)};
79
80 // Construct an extended pdf, which measures the event count N **on the full range**.
81 // If the actual domain of x that is fitted is identical to FULL, this has no affect.
82 //
83 // If the fitted domain is a subset of FULL, though, the expected event count is divided by
84 // \f[
85 // \mathrm{frac} = \frac{
86 // \int_{\mathrm{Fit range}} \mathrm{model}(x) \; \mathrm{d}x }{
87 // \int_{\mathrm{Full range}} \mathrm{model}(x) \; \mathrm{d}x }.
88 // \f]
89 // N will therefore return the count extrapolated to the full range instead of the fit range.
90 //
91 // **Note**: When using a RooAddPdf for extending the likelihood, the same effect can be achieved with
93
94 RooRealVar N("N", "Extended term", 0, 20000);
95 RooExtendPdf extmodel("extmodel", "Extended model", model, N, "FULL");
96
97
98 // It can be instructive to fit the above model to either the LEFT or RIGHT range. N should approximately converge to the expected number
99 // of events in the full range. One may try to leave out "FULL" in the constructor, or the interpretation of N changes.
100 extmodel.fitTo(*data, Range("LEFT"), PrintLevel(-1));
101 N.Print();
102
103
104 // If we now do a simultaneous fit to the extended model, instead of the original model, the LEFT and RIGHT range will each correct their local
105 // N such that it refers to the FULL range.
106 //
107 // This joint fit of the extmodel will include (w.r.t. the plain model fit) a product of extended terms
108 // \f[
109 // L_\mathrm{ext} = L
110 // \cdot \mathrm{Poisson} \left( N_\mathrm{obs}^\mathrm{LEFT} | N_\mathrm{exp} / \mathrm{frac LEFT} \right)
111 // \cdot \mathrm{Poisson} \left( N_\mathrm{obs}^\mathrm{RIGHT} | N_\mathrm{exp} / \mathrm{frac RIGHT} \right)
112 // \f]
113
114
115 TCanvas* c = new TCanvas("c", "c", 2100, 700);
116 c->Divide(3);
117 c->cd(1);
118
119 std::unique_ptr<RooFitResult> r{model.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())};
120 r->Print();
121
122 RooPlot* frame = x.frame();
123 data->plotOn(frame);
124 model.plotOn(frame, VisualizeError(*r));
125 model.plotOn(frame);
126 model.paramOn(frame, Label("Non-extended fit"));
127 frame->Draw();
128
129 c->cd(2);
130
131 std::unique_ptr<RooFitResult> r2{extmodel.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())};
132 r2->Print();
133 RooPlot* frame2 = x.frame();
134 data->plotOn(frame2);
135 extmodel.plotOn(frame2);
136 extmodel.plotOn(frame2, VisualizeError(*r2));
137 extmodel.paramOn(frame2, Label("Extended fit"), Layout(0.4,0.95));
138 frame2->Draw();
139
140 // PART 2: Extending with RooAddPdf
141 // --------------------------------
142 //
143 // Now we repeat the above exercise, but instead of explicitly adding an extended term to a single shape pdf (RooExponential),
144 // we assume that we already have an extended likelihood model in the form of a RooAddPdf constructed in the form Nsig * sigPdf + Nbkg * bkgPdf.
145 //
146 // We add a Gaussian to the previously defined exponential background.
147 // The signal shape parameters are chosen constant, since the signal region is entirely in the blinded region, i.e., the fit has no sensitivity.
148
149 c->cd(3);
150
151 RooRealVar Nsig("Nsig", "Number of signal events", 1000, 0, 2000);
152 RooRealVar Nbkg("Nbkg", "Number of background events", 10000, 0, 20000);
153
154 RooRealVar mean("mean", "Mean of signal model", 40.);
155 RooRealVar width("width", "Width of signal model", 5.);
156 RooGaussian sig("sig", "Signal model", x, mean, width);
157
158 RooAddPdf modelsum("modelsum", "NSig*signal + NBkg*background", {sig, model}, {Nsig, Nbkg});
159
160 // This model will automatically insert the correction factor for the reinterpretation of Nsig and Nnbkg in the full ranges.
161 //
162 // When this happens, it reports this with lines like the following:
163 // 
164 // [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_LEFT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables
165 // [#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_RIGHT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables
166 // 
167
168 std::unique_ptr<RooFitResult> r3{modelsum.fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1), Save())};
169 r3->Print();
170
171 RooPlot* frame3 = x.frame();
172 data->plotOn(frame3);
173 modelsum.plotOn(frame3);
174 modelsum.plotOn(frame3, VisualizeError(*r3));
175 modelsum.paramOn(frame3, Label("S+B fit with RooAddPdf"), Layout(0.3,0.95));
176 frame3->Draw();
177
178 c->Draw();
179}
#define c(i)
Definition RSha256.hxx:101
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t width
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
Efficient implementation of a sum of PDFs of the form.
Exponential PDF.
RooExtendPdf is a wrapper around an existing PDF that adds a parameteric extended likelihood term to ...
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:237
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:649
Variable that can be changed from the outside.
Definition RooRealVar.h:37
The Canvas class.
Definition TCanvas.h:23
RooCmdArg Layout(double xmin, double xmax=0.99, double ymin=0.95)
RooCmdArg Save(bool flag=true)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg VisualizeError(const RooDataSet &paramData, double Z=1)
Double_t x[n]
Definition legend1.C:17
const Rcpp::internal::NamedPlaceHolder & Label
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Ta Range(0, 0, 1, 1)