Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs301_splot.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// SPlot tutorial
5///
6/// This tutorial shows an example of using SPlot to unfold two distributions.
7/// The physics context for the example is that we want to know
8/// the isolation distribution for real electrons from Z events
9/// and fake electrons from QCD. Isolation is our 'control' variable.
10/// To unfold them, we need a model for an uncorrelated variable that
11/// discriminates between Z and QCD. To do this, we use the invariant
12/// mass of two electrons. We model the Z with a Gaussian and the QCD
13/// with a falling exponential.
14///
15/// Note, since we don't have real data in this tutorial, we need to generate
16/// toy data. To do that we need a model for the isolation variable for
17/// both Z and QCD. This is only used to generate the toy data, and would
18/// not be needed if we had real data.
19///
20/// \macro_image
21/// \macro_code
22/// \macro_output
23///
24/// \author Kyle Cranmer
25
26#include "RooRealVar.h"
27#include "RooStats/SPlot.h"
28#include "RooDataSet.h"
29#include "RooRealVar.h"
30#include "RooGaussian.h"
31#include "RooExponential.h"
32#include "RooChebychev.h"
33#include "RooAddPdf.h"
34#include "RooProdPdf.h"
35#include "RooAddition.h"
36#include "RooProduct.h"
37#include "TCanvas.h"
38#include "RooAbsPdf.h"
39#include "RooFit.h"
40#include "RooFitResult.h"
41#include "RooWorkspace.h"
42#include "RooConstVar.h"
43#include <iomanip>
44
45// use this order for safety on library loading
46using namespace RooFit;
47using namespace RooStats;
48
49// see below for implementation
50void AddModel(RooWorkspace *);
51void AddData(RooWorkspace *);
52void DoSPlot(RooWorkspace *);
53void MakePlots(RooWorkspace *);
54
55void rs301_splot()
56{
57
58 // Create a new workspace to manage the project.
59 RooWorkspace *wspace = new RooWorkspace("myWS");
60
61 // add the signal and background models to the workspace.
62 // Inside this function you will find a description of our model.
63 AddModel(wspace);
64
65 // add some toy data to the workspace
66 AddData(wspace);
67
68 // inspect the workspace if you wish
69 // wspace->Print();
70
71 // do sPlot.
72 // This will make a new dataset with sWeights added for every event.
73 DoSPlot(wspace);
74
75 // Make some plots showing the discriminating variable and
76 // the control variable after unfolding.
77 MakePlots(wspace);
78
79 // cleanup
80 delete wspace;
81}
82
83//____________________________________
84void AddModel(RooWorkspace *ws)
85{
86
87 // Make models for signal (Higgs) and background (Z+jets and QCD)
88 // In real life, this part requires an intelligent modeling
89 // of signal and background -- this is only an example.
90
91 // set range of observable
92 Double_t lowRange = 0., highRange = 200.;
93
94 // make a RooRealVar for the observables
95 RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange, "GeV");
96 RooRealVar isolation("isolation", "isolation", 0., 20., "GeV");
97
98 // --------------------------------------
99 // make 2-d model for Z including the invariant mass
100 // distribution and an isolation distribution which we want to
101 // unfold from QCD.
102 std::cout << "make z model" << std::endl;
103 // mass model for Z
104 RooRealVar mZ("mZ", "Z Mass", 91.2, lowRange, highRange);
105 RooRealVar sigmaZ("sigmaZ", "Width of Gaussian", 2, 0, 10, "GeV");
106 RooGaussian mZModel("mZModel", "Z+jets Model", invMass, mZ, sigmaZ);
107 // we know Z mass
108 mZ.setConstant();
109 // we leave the width of the Z free during the fit in this example.
110
111 // isolation model for Z. Only used to generate toy MC.
112 // the exponential is of the form exp(c*x). If we want
113 // the isolation to decay an e-fold every R GeV, we use
114 // c = -1/R.
115 RooConstVar zIsolDecayConst("zIsolDecayConst", "z isolation decay constant", -1);
116 RooExponential zIsolationModel("zIsolationModel", "z isolation model", isolation, zIsolDecayConst);
117
118 // make the combined Z model
119 RooProdPdf zModel("zModel", "2-d model for Z", RooArgSet(mZModel, zIsolationModel));
120
121 // --------------------------------------
122 // make QCD model
123
124 std::cout << "make qcd model" << std::endl;
125 // mass model for QCD.
126 // the exponential is of the form exp(c*x). If we want
127 // the mass to decay an e-fold every R GeV, we use
128 // c = -1/R.
129 // We can leave this parameter free during the fit.
130 RooRealVar qcdMassDecayConst("qcdMassDecayConst", "Decay const for QCD mass spectrum", -0.01, -100, 100, "1/GeV");
131 RooExponential qcdMassModel("qcdMassModel", "qcd Mass Model", invMass, qcdMassDecayConst);
132
133 // isolation model for QCD. Only used to generate toy MC
134 // the exponential is of the form exp(c*x). If we want
135 // the isolation to decay an e-fold every R GeV, we use
136 // c = -1/R.
137 RooConstVar qcdIsolDecayConst("qcdIsolDecayConst", "Et resolution constant", -.1);
138 RooExponential qcdIsolationModel("qcdIsolationModel", "QCD isolation model", isolation, qcdIsolDecayConst);
139
140 // make the 2-d model
141 RooProdPdf qcdModel("qcdModel", "2-d model for QCD", RooArgSet(qcdMassModel, qcdIsolationModel));
142
143 // --------------------------------------
144 // combined model
145
146 // These variables represent the number of Z or QCD events
147 // They will be fitted.
148 RooRealVar zYield("zYield", "fitted yield for Z", 50, 0., 1000);
149 RooRealVar qcdYield("qcdYield", "fitted yield for QCD", 100, 0., 1000);
150
151 // now make the combined model
152 std::cout << "make full model" << std::endl;
153 RooAddPdf model("model", "z+qcd background models", RooArgList(zModel, qcdModel), RooArgList(zYield, qcdYield));
154
155 // interesting for debugging and visualizing the model
156 model.graphVizTree("fullModel.dot");
157
158 std::cout << "import model" << std::endl;
159
160 ws->import(model);
161}
162
163//____________________________________
164void AddData(RooWorkspace *ws)
165{
166 // Add a toy dataset
167
168 // how many events do we want?
169 Int_t nEvents = 1000;
170
171 // get what we need out of the workspace to make toy data
172 RooAbsPdf *model = ws->pdf("model");
173 RooRealVar *invMass = ws->var("invMass");
174 RooRealVar *isolation = ws->var("isolation");
175
176 // make the toy data
177 std::cout << "make data set and import to workspace" << std::endl;
178 RooDataSet *data = model->generate(RooArgSet(*invMass, *isolation), nEvents);
179
180 // import data into workspace
181 ws->import(*data, Rename("data"));
182}
183
184//____________________________________
185void DoSPlot(RooWorkspace *ws)
186{
187 std::cout << "Calculate sWeights" << std::endl;
188
189 // get what we need out of the workspace to do the fit
190 RooAbsPdf *model = ws->pdf("model");
191 RooRealVar *zYield = ws->var("zYield");
192 RooRealVar *qcdYield = ws->var("qcdYield");
193 RooDataSet *data = (RooDataSet *)ws->data("data");
194
195 // fit the model to the data.
196 model->fitTo(*data, Extended());
197
198 // The sPlot technique requires that we fix the parameters
199 // of the model that are not yields after doing the fit.
200 //
201 // This *could* be done with the lines below, however this is taken care of
202 // by the RooStats::SPlot constructor (or more precisely the AddSWeight
203 // method).
204 //
205 // RooRealVar* sigmaZ = ws->var("sigmaZ");
206 // RooRealVar* qcdMassDecayConst = ws->var("qcdMassDecayConst");
207 // sigmaZ->setConstant();
208 // qcdMassDecayConst->setConstant();
209
211
212 std::cout << "\n\n------------------------------------------\nThe dataset before creating sWeights:\n";
213 data->Print();
214
216
217 // Now we use the SPlot class to add SWeights to our data set
218 // based on our model and our yield variables
219 RooStats::SPlot *sData = new RooStats::SPlot("sData", "An SPlot", *data, model, RooArgList(*zYield, *qcdYield));
220
221 std::cout << "\n\nThe dataset after creating sWeights:\n";
222 data->Print();
223
224 // Check that our weights have the desired properties
225
226 std::cout << "\n\n------------------------------------------\n\nCheck SWeights:" << std::endl;
227
228 std::cout << std::endl
229 << "Yield of Z is\t" << zYield->getVal() << ". From sWeights it is "
230 << sData->GetYieldFromSWeight("zYield") << std::endl;
231
232 std::cout << "Yield of QCD is\t" << qcdYield->getVal() << ". From sWeights it is "
233 << sData->GetYieldFromSWeight("qcdYield") << std::endl
234 << std::endl;
235
236 for (Int_t i = 0; i < 10; i++) {
237 std::cout << "z Weight for event " << i << std::right << std::setw(12) << sData->GetSWeight(i, "zYield") << " qcd Weight"
238 << std::setw(12) << sData->GetSWeight(i, "qcdYield") << " Total Weight" << std::setw(12) << sData->GetSumOfEventSWeight(i)
239 << std::endl;
240 }
241
242 std::cout << std::endl;
243
244 // import this new dataset with sWeights
245 std::cout << "import new dataset with sWeights" << std::endl;
246 ws->import(*data, Rename("dataWithSWeights"));
247
249}
250
251void MakePlots(RooWorkspace *ws)
252{
253
254 // Here we make plots of the discriminating variable (invMass) after the fit
255 // and of the control variable (isolation) after unfolding with sPlot.
256 std::cout << "make plots" << std::endl;
257
258 // make our canvas
259 TCanvas *cdata = new TCanvas("sPlot", "sPlot demo", 400, 600);
260 cdata->Divide(1, 3);
261
262 // get what we need out of the workspace
263 RooAbsPdf *model = ws->pdf("model");
264 RooAbsPdf *zModel = ws->pdf("zModel");
265 RooAbsPdf *qcdModel = ws->pdf("qcdModel");
266
267 RooRealVar *isolation = ws->var("isolation");
268 RooRealVar *invMass = ws->var("invMass");
269
270 // note, we get the dataset with sWeights
271 RooDataSet *data = (RooDataSet *)ws->data("dataWithSWeights");
272
273 // this shouldn't be necessary, need to fix something with workspace
274 // do this to set parameters back to their fitted values.
275// model->fitTo(*data, Extended());
276
277 // plot invMass for data with full model and individual components overlaid
278 // TCanvas* cdata = new TCanvas();
279 cdata->cd(1);
280 RooPlot *frame = invMass->frame();
281 data->plotOn(frame);
282 model->plotOn(frame, Name("FullModel"));
283 model->plotOn(frame, Components(*zModel), LineStyle(kDashed), LineColor(kRed), Name("ZModel"));
284 model->plotOn(frame, Components(*qcdModel), LineStyle(kDashed), LineColor(kGreen), Name("QCDModel"));
285
286 TLegend leg(0.11, 0.5, 0.5, 0.8);
287 leg.AddEntry(frame->findObject("FullModel"), "Full model", "L");
288 leg.AddEntry(frame->findObject("ZModel"), "Z model", "L");
289 leg.AddEntry(frame->findObject("QCDModel"), "QCD model", "L");
290 leg.SetBorderSize(0);
291 leg.SetFillStyle(0);
292
293 frame->SetTitle("Fit of model to discriminating variable");
294 frame->Draw();
295 leg.DrawClone();
296
297 // Now use the sWeights to show isolation distribution for Z and QCD.
298 // The SPlot class can make this easier, but here we demonstrate in more
299 // detail how the sWeights are used. The SPlot class should make this
300 // very easy and needs some more development.
301
302 // Plot isolation for Z component.
303 // Do this by plotting all events weighted by the sWeight for the Z component.
304 // The SPlot class adds a new variable that has the name of the corresponding
305 // yield + "_sw".
306 cdata->cd(2);
307
308 // create weighted data set
309 RooDataSet *dataw_z = new RooDataSet(data->GetName(), data->GetTitle(), data, *data->get(), 0, "zYield_sw");
310
311 RooPlot *frame2 = isolation->frame();
312 // Since the data are weighted, we use SumW2 to compute the errors.
313 dataw_z->plotOn(frame2, DataError(RooAbsData::SumW2));
314
315 frame2->SetTitle("Isolation distribution with s weights to project out Z");
316 frame2->Draw();
317
318 // Plot isolation for QCD component.
319 // Eg. plot all events weighted by the sWeight for the QCD component.
320 // The SPlot class adds a new variable that has the name of the corresponding
321 // yield + "_sw".
322 cdata->cd(3);
323 RooDataSet *dataw_qcd = new RooDataSet(data->GetName(), data->GetTitle(), data, *data->get(), 0, "qcdYield_sw");
324 RooPlot *frame3 = isolation->frame();
325 dataw_qcd->plotOn(frame3, DataError(RooAbsData::SumW2));
326
327 frame3->SetTitle("Isolation distribution with s weights to project out QCD");
328 frame3->Draw();
329
330 // cdata->SaveAs("SPlot.gif");
331}
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kDashed
Definition TAttLine.h:48
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition RooAbsData.h:193
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:58
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:121
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:91
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
Definition RooAddPdf.h:32
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
RooConstVar represent a constant real-valued object.
Definition RooConstVar.h:26
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:33
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Exponential PDF.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
void setSilentMode(Bool_t flag)
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:44
void SetTitle(const char *name)
Set the title of the RooPlot to 'title'.
Definition RooPlot.cxx:1242
TObject * findObject(const char *name, const TClass *clas=0) const
Find the named object in our list of items and return a pointer to it.
Definition RooPlot.cxx:982
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:691
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:37
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
A class to calculate "sWeights" used to create an "sPlot".
Definition SPlot.h:32
Double_t GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition SPlot.cxx:331
Double_t GetSWeight(Int_t numEvent, const char *sVariable) const
Retrieve an s weight.
Definition SPlot.cxx:251
Double_t GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition SPlot.cxx:300
The RooWorkspace is a persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:708
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
virtual const char * GetTitle() const
Returns title of object.
Definition TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1177
RooCmdArg Rename(const char *suffix)
RooCmdArg Extended(Bool_t flag=kTRUE)
RooCmdArg DataError(Int_t)
RooCmdArg Components(const RooArgSet &compSet)
RooCmdArg LineColor(Color_t color)
RooCmdArg LineStyle(Style_t style)
leg
Definition legend1.C:34
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition Asimov.h:19
const char * Name
Definition TXMLSetup.cxx:67
void ws()
Definition ws.C:66