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