Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs102_hypotestwithshapes.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// A typical search for a new particle by studying an invariant mass distribution
5///
6/// The macro creates a simple signal model and two background models,
7/// which are added to a RooWorkspace.
8/// The macro creates a toy dataset, and then uses a RooStats
9/// ProfileLikleihoodCalculator to do a hypothesis test of the
10/// background-only and signal+background hypotheses.
11/// In this example, shape uncertainties are not taken into account, but
12/// normalization uncertainties are.
13///
14/// \macro_image
15/// \macro_output
16/// \macro_code
17///
18/// \author Kyle Cranmer
19
20#include "RooDataSet.h"
21#include "RooRealVar.h"
22#include "RooGaussian.h"
23#include "RooAddPdf.h"
24#include "RooProdPdf.h"
25#include "RooAddition.h"
26#include "RooProduct.h"
27#include "TCanvas.h"
28#include "RooChebychev.h"
29#include "RooAbsPdf.h"
30#include "RooFitResult.h"
31#include "RooPlot.h"
32#include "RooAbsArg.h"
33#include "RooWorkspace.h"
36#include <string>
37
38// use this order for safety on library loading
39using namespace RooFit;
40using namespace RooStats;
41
42// see below for implementation
43void AddModel(RooWorkspace *);
44void AddData(RooWorkspace *);
45void DoHypothesisTest(RooWorkspace *);
46void MakePlots(RooWorkspace *);
47
48//____________________________________
49void rs102_hypotestwithshapes()
50{
51
52 // The main macro.
53
54 // Create a workspace to manage the project.
55 RooWorkspace *wspace = new RooWorkspace("myWS");
56
57 // add the signal and background models to the workspace
58 AddModel(wspace);
59
60 // add some toy data to the workspace
61 AddData(wspace);
62
63 // inspect the workspace if you wish
64 // wspace->Print();
65
66 // do the hypothesis test
67 DoHypothesisTest(wspace);
68
69 // make some plots
70 MakePlots(wspace);
71
72 // cleanup
73 delete wspace;
74}
75
76//____________________________________
77void AddModel(RooWorkspace *wks)
78{
79
80 // Make models for signal (Higgs) and background (Z+jets and QCD)
81 // In real life, this part requires an intelligent modeling
82 // of signal and background -- this is only an example.
83
84 // set range of observable
85 Double_t lowRange = 60, highRange = 200;
86
87 // make a RooRealVar for the observable
88 RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange, "GeV");
89
90 // --------------------------------------
91 // make a simple signal model.
92 RooRealVar mH("mH", "Higgs Mass", 130, 90, 160);
93 RooRealVar sigma1("sigma1", "Width of Gaussian", 12., 2, 100);
94 RooGaussian sigModel("sigModel", "Signal Model", invMass, mH, sigma1);
95 // we will test this specific mass point for the signal
96 mH.setConstant();
97 // and we assume we know the mass resolution
98 sigma1.setConstant();
99
100 // --------------------------------------
101 // make zjj model. Just like signal model
102 RooRealVar mZ("mZ", "Z Mass", 91.2, 0, 100);
103 RooRealVar sigma1_z("sigma1_z", "Width of Gaussian", 10., 6, 100);
104 RooGaussian zjjModel("zjjModel", "Z+jets Model", invMass, mZ, sigma1_z);
105 // we know Z mass
106 mZ.setConstant();
107 // assume we know resolution too
108 sigma1_z.setConstant();
109
110 // --------------------------------------
111 // make QCD model
112 RooRealVar a0("a0", "a0", 0.26, -1, 1);
113 RooRealVar a1("a1", "a1", -0.17596, -1, 1);
114 RooRealVar a2("a2", "a2", 0.018437, -1, 1);
115 RooRealVar a3("a3", "a3", 0.02, -1, 1);
116 RooChebychev qcdModel("qcdModel", "A Polynomial for QCD", invMass, RooArgList(a0, a1, a2));
117
118 // let's assume this shape is known, but the normalization is not
119 a0.setConstant();
120 a1.setConstant();
121 a2.setConstant();
122
123 // --------------------------------------
124 // combined model
125
126 // Setting the fraction of Zjj to be 40% for initial guess.
127 RooRealVar fzjj("fzjj", "fraction of zjj background events", .4, 0., 1);
128
129 // Set the expected fraction of signal to 20%.
130 RooRealVar fsigExpected("fsigExpected", "expected fraction of signal events", .2, 0., 1);
131 fsigExpected.setConstant(); // use mu as main parameter, so fix this.
132
133 // Introduce mu: the signal strength in units of the expectation.
134 // eg. mu = 1 is the SM, mu = 0 is no signal, mu=2 is 2x the SM
135 RooRealVar mu("mu", "signal strength in units of SM expectation", 1, 0., 2);
136
137 // Introduce ratio of signal efficiency to nominal signal efficiency.
138 // This is useful if you want to do limits on cross section.
139 RooRealVar ratioSigEff("ratioSigEff", "ratio of signal efficiency to nominal signal efficiency", 1., 0., 2);
140 ratioSigEff.setConstant(kTRUE);
141
142 // finally the signal fraction is the product of the terms above.
143 RooProduct fsig("fsig", "fraction of signal events", RooArgSet(mu, ratioSigEff, fsigExpected));
144
145 // full model
146 RooAddPdf model("model", "sig+zjj+qcd background shapes", RooArgList(sigModel, zjjModel, qcdModel),
147 RooArgList(fsig, fzjj));
148
149 // interesting for debugging and visualizing the model
150 // model.printCompactTree("","fullModel.txt");
151 // model.graphVizTree("fullModel.dot");
152
153 wks->import(model);
154}
155
156//____________________________________
157void AddData(RooWorkspace *wks)
158{
159 // Add a toy dataset
160
161 Int_t nEvents = 150;
162 RooAbsPdf *model = wks->pdf("model");
163 RooRealVar *invMass = wks->var("invMass");
164
165 std::unique_ptr<RooDataSet> data{model->generate(*invMass, nEvents)};
166
167 wks->import(*data, Rename("data"));
168}
169
170//____________________________________
171void DoHypothesisTest(RooWorkspace *wks)
172{
173
174 // Use a RooStats ProfileLikleihoodCalculator to do the hypothesis test.
175 ModelConfig model;
176 model.SetWorkspace(*wks);
177 model.SetPdf("model");
178
179 // plc.SetData("data");
180
182 plc.SetData(*(wks->data("data")));
183
184 // here we explicitly set the value of the parameters for the null.
185 // We want no signal contribution, eg. mu = 0
186 RooRealVar *mu = wks->var("mu");
187 // RooArgSet* nullParams = new RooArgSet("nullParams");
188 // nullParams->addClone(*mu);
189 RooArgSet poi(*mu);
190 RooArgSet *nullParams = (RooArgSet *)poi.snapshot();
191 nullParams->setRealValue("mu", 0);
192
193 // plc.SetNullParameters(*nullParams);
194 plc.SetModel(model);
195 // NOTE: using snapshot will import nullparams
196 // in the WS and merge with existing "mu"
197 // model.SetSnapshot(*nullParams);
198
199 // use instead setNuisanceParameters
200 plc.SetNullParameters(*nullParams);
201
202 // We get a HypoTestResult out of the calculator, and we can query it.
203 HypoTestResult *htr = plc.GetHypoTest();
204 cout << "-------------------------------------------------" << endl;
205 cout << "The p-value for the null is " << htr->NullPValue() << endl;
206 cout << "Corresponding to a significance of " << htr->Significance() << endl;
207 cout << "-------------------------------------------------\n\n" << endl;
208}
209
210//____________________________________
211void MakePlots(RooWorkspace *wks)
212{
213
214 // Make plots of the data and the best fit model in two cases:
215 // first the signal+background case
216 // second the background-only case.
217
218 // get some things out of workspace
219 RooAbsPdf *model = wks->pdf("model");
220 RooAbsPdf *sigModel = wks->pdf("sigModel");
221 RooAbsPdf *zjjModel = wks->pdf("zjjModel");
222 RooAbsPdf *qcdModel = wks->pdf("qcdModel");
223
224 RooRealVar *mu = wks->var("mu");
225 RooRealVar *invMass = wks->var("invMass");
226 RooAbsData *data = wks->data("data");
227
228 // --------------------------------------
229 // Make plots for the Alternate hypothesis, eg. let mu float
230
231 mu->setConstant(kFALSE);
232
233 model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE), PrintLevel(-1));
234
235 // plot sig candidates, full model, and individual components
236 new TCanvas();
237 RooPlot *frame = invMass->frame();
238 data->plotOn(frame);
239 model->plotOn(frame);
240 model->plotOn(frame, Components(*sigModel), LineStyle(kDashed), LineColor(kRed));
241 model->plotOn(frame, Components(*zjjModel), LineStyle(kDashed), LineColor(kBlack));
242 model->plotOn(frame, Components(*qcdModel), LineStyle(kDashed), LineColor(kGreen));
243
244 frame->SetTitle("An example fit to the signal + background model");
245 frame->Draw();
246 // cdata->SaveAs("alternateFit.gif");
247
248 // --------------------------------------
249 // Do Fit to the Null hypothesis. Eg. fix mu=0
250
251 mu->setVal(0); // set signal fraction to 0
252 mu->setConstant(kTRUE); // set constant
253
254 model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE), PrintLevel(-1));
255
256 // plot signal candidates with background model and components
257 new TCanvas();
258 RooPlot *xframe2 = invMass->frame();
259 data->plotOn(xframe2, DataError(RooAbsData::SumW2));
260 model->plotOn(xframe2);
261 model->plotOn(xframe2, Components(*zjjModel), LineStyle(kDashed), LineColor(kBlack));
262 model->plotOn(xframe2, Components(*qcdModel), LineStyle(kDashed), LineColor(kGreen));
263
264 xframe2->SetTitle("An example fit to the background-only model");
265 xframe2->Draw();
266 // cbkgonly->SaveAs("nullFit.gif");
267}
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
@ kRed
Definition Rtypes.h:66
@ kBlack
Definition Rtypes.h:65
@ kGreen
Definition Rtypes.h:66
@ kDashed
Definition TAttLine.h:48
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
bool setRealValue(const char *name, double newVal=0.0, bool verbose=false)
Set value of a RooAbsRealLValue stored in set with given name to newVal No error messages are printed...
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:124
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:57
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
void setConstant(bool value=true)
Efficient implementation of a sum of PDFs of the form.
Definition RooAddPdf.h:33
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:154
Chebychev polynomial p.d.f.
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
void SetTitle(const char *name) override
Set the title of the RooPlot to 'title'.
Definition RooPlot.cxx:1203
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:597
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
virtual void SetNullParameters(const RooArgSet &set)
set parameter values for the null if using a common PDF
void SetModel(const ModelConfig &model) override
set the model (in this case can set only the parameters for the null hypothesis)
void SetData(RooAbsData &data) override
Set the DataSet, add to the workspace if not already there.
HypoTestResult is a base class for results from hypothesis tests.
virtual double Significance() const
familiar name for the Null p-value in terms of 1-sided Gaussian significance
virtual double NullPValue() const
Return p-value for null hypothesis.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
virtual void SetWorkspace(RooWorkspace &ws)
Definition ModelConfig.h:71
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the workspace if not already there.
Definition ModelConfig.h:88
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
HypoTestResult * GetHypoTest() const override
Return the hypothesis test result obtained from the likelihood ratio of the maximum likelihood value ...
Persistable container for RooFit projects.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
RooAbsData * data(RooStringView name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
The Canvas class.
Definition TCanvas.h:23
RooCmdArg Rename(const char *suffix)
RooCmdArg Hesse(bool flag=true)
RooCmdArg Save(bool flag=true)
RooCmdArg Minos(bool flag=true)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg DataError(Int_t)
RooCmdArg Components(Args_t &&... argsOrArgSet)
RooCmdArg LineColor(Color_t color)
RooCmdArg LineStyle(Style_t style)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:64
Namespace for the RooStats classes.
Definition CodegenImpl.h:58