Logo ROOT   master
Reference Guide
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
46 using namespace RooFit;
47 using namespace RooStats;
48 
49 // see below for implementation
50 void AddModel(RooWorkspace *);
51 void AddData(RooWorkspace *);
52 void DoSPlot(RooWorkspace *);
53 void MakePlots(RooWorkspace *);
54 
55 void 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 //____________________________________
84 void 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 //____________________________________
164 void 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 //____________________________________
185 void 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 
251 void 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 }
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
Definition: RooAbsPdf.h:55
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
Definition: RooAddPdf.h:29
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
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
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:546
void ws()
Definition: ws.C:66
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooProdPdf is an efficient implementation of a product of PDFs of the form .
Definition: RooProdPdf.h:31
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:697
RooCmdArg DataError(Int_t)
void SetTitle(const char *name)
Set the title of the RooPlot to &#39;title&#39;.
Definition: RooPlot.cxx:1257
Definition: Rtypes.h:64
static RooMsgService & instance()
Return reference to singleton instance.
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:118
Double_t GetYieldFromSWeight(const char *sVariable) const
Sum the SWeights for a particular species over all events.
Definition: SPlot.cxx:280
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
Double_t GetSWeight(Int_t numEvent, const char *sVariable) const
Retrieve an s weight.
Definition: SPlot.cxx:200
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
Exponential PDF.
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
RooCmdArg Name(const char *name)
RooCmdArg Rename(const char *suffix)
void setGlobalKillBelow(RooFit::MsgLevel level)
void setSilentMode(Bool_t flag)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
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, but no plot contents.
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The Canvas class.
Definition: TCanvas.h:27
Namespace for the RooStats classes.
Definition: Asimov.h:19
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
RooCmdArg Extended(Bool_t flag=kTRUE)
leg
Definition: legend1.C:34
RooCmdArg LineColor(Color_t color)
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1164
RooCmdArg LineStyle(Style_t style)
RooCmdArg Components(const RooArgSet &compSet)
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.
Definition: RooAbsPdf.cxx:1254
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Double_t GetSumOfEventSWeight(Int_t numEvent) const
Sum the SWeights for a particular event.
Definition: SPlot.cxx:249
This class calculates sWeights used to create an sPlot.
Definition: SPlot.h:32
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
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:1002
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:711
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48