Logo ROOT  
Reference Guide
RooStatsUtils.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer 28/07/2008
3
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#include "Rtypes.h"
13
14#if !defined(R__ALPHA) && !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
16#endif
17
18#include "TTree.h"
19#include "TBranch.h"
20
21#include "RooUniform.h"
22#include "RooProdPdf.h"
23#include "RooExtendPdf.h"
24#include "RooSimultaneous.h"
27#include <typeinfo>
28
29using namespace std;
30
31namespace RooStats {
32
34 static RooStatsConfig theConfig;
35 return theConfig;
36 }
37
39 // Asimov significance
40 // formula [10] and [20] from https://www.pp.rhul.ac.uk/~cowan/stat/notes/medsigNote.pdf
41 // case we have a sigma_b
42 double sb2 = sigma_b*sigma_b;
43 // formula below has a large error when sigma_b becomes zero
44 // better to use the approximation for sigma_b=0 for very small values
45 double r = sb2/b;
46 if (r > 1.E-12) {
47 double bpsb2 = b + sb2;
48 double b2 = b*b;
49 double spb = s+b;
50 double za2 = 2.*( (spb)* std::log( ( spb)*(bpsb2)/(b2+ spb*sb2) ) -
51 (b2/sb2) * std::log(1. + ( sb2 * s)/(b * bpsb2) ) );
52 return sqrt(za2);
53
54 }
55 // case when the background (b) is known
56 double za2 = 2.*( (s+b) * std::log(1. + s/b) -s );
57 return std::sqrt(za2);
58 }
59
60 /// Use an offset in NLL calculations.
61 void UseNLLOffset(bool on) {
63 }
64
65 /// Test of RooStats should by default offset NLL calculations.
66 bool IsNLLOffset() {
68 }
69
70 void FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
71 // utility function to factorize constraint terms from a pdf
72 // (from G. Petrucciani)
73 const std::type_info & id = typeid(pdf);
74 if (id == typeid(RooProdPdf)) {
75 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
76 RooArgList list(prod->pdfList());
77 for (int i = 0, n = list.getSize(); i < n; ++i) {
78 RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
79 FactorizePdf(observables, *pdfi, obsTerms, constraints);
80 }
81 } else if (id == typeid(RooExtendPdf)) {
82 TIterator *iter = pdf.serverIterator();
83 // extract underlying pdf which is extended; first server is the pdf; second server is the number of events variable
84 RooAbsPdf *updf = dynamic_cast<RooAbsPdf *>(iter->Next());
85 assert(updf != 0);
86 delete iter;
87 FactorizePdf(observables, *updf, obsTerms, constraints);
88 } else if (id == typeid(RooSimultaneous)) { //|| id == typeid(RooSimultaneousOpt)) {
89 RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf);
90 assert(sim != 0);
92 for (int ic = 0, nc = cat->numBins((const char *)0); ic < nc; ++ic) {
93 cat->setBin(ic);
94 RooAbsPdf* catPdf = sim->getPdf(cat->getCurrentLabel());
95 // it is possible that a pdf is not defined for every category
96 if (catPdf != 0) FactorizePdf(observables, *catPdf, obsTerms, constraints);
97 }
98 delete cat;
99 } else if (pdf.dependsOn(observables)) {
100 if (!obsTerms.contains(pdf)) obsTerms.add(pdf);
101 } else {
102 if (!constraints.contains(pdf)) constraints.add(pdf);
103 }
104 }
105
106
107 void FactorizePdf(RooStats::ModelConfig &model, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
108 // utility function to factorize constraint terms from a pdf
109 // (from G. Petrucciani)
110 if (!model.GetObservables() ) {
111 oocoutE((TObject*)0,InputArguments) << "RooStatsUtils::FactorizePdf - invalid input model: missing observables" << endl;
112 return;
113 }
114 return FactorizePdf(*model.GetObservables(), pdf, obsTerms, constraints);
115 }
116
117
118 RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name) {
119 // make a nuisance pdf by factorizing out all constraint terms in a common pdf
120 RooArgList obsTerms, constraints;
121 FactorizePdf(observables, pdf, obsTerms, constraints);
122 if(constraints.getSize() == 0) {
123 oocoutW((TObject *)0, Eval) << "RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model" << endl;
124 return 0;
125 }
126 return new RooProdPdf(name,"", constraints);
127 }
128
130 // make a nuisance pdf by factorizing out all constraint terms in a common pdf
131 if (!model.GetPdf() || !model.GetObservables() ) {
132 oocoutE((TObject*)0, InputArguments) << "RooStatsUtils::MakeNuisancePdf - invalid input model: missing pdf and/or observables" << endl;
133 return 0;
134 }
135 return MakeNuisancePdf(*model.GetPdf(), *model.GetObservables(), name);
136 }
137
138 RooAbsPdf * StripConstraints(RooAbsPdf &pdf, const RooArgSet &observables) {
139 const std::type_info & id = typeid(pdf);
140
141 if (id == typeid(RooProdPdf)) {
142
143 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
144 RooArgList list(prod->pdfList()); RooArgList newList;
145
146 for (int i = 0, n = list.getSize(); i < n; ++i) {
147 RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
148 RooAbsPdf *newPdfi = StripConstraints(*pdfi, observables);
149 if(newPdfi != NULL) newList.add(*newPdfi);
150 }
151
152 if(newList.getSize() == 0) return NULL; // only constraints in product
153 // return single component (no longer a product)
154 else if(newList.getSize() == 1) return dynamic_cast<RooAbsPdf *>(newList.at(0)->clone(TString::Format("%s_unconstrained",
155 newList.at(0)->GetName())));
156 else return new RooProdPdf(TString::Format("%s_unconstrained", prod->GetName()).Data(),
157 TString::Format("%s without constraints", prod->GetTitle()).Data(), newList);
158
159 } else if (id == typeid(RooExtendPdf)) {
160
161 TIterator *iter = pdf.serverIterator();
162 // extract underlying pdf which is extended; first server is the pdf; second server is the number of events variable
163 RooAbsPdf *uPdf = dynamic_cast<RooAbsPdf *>(iter->Next());
164 RooAbsReal *extended_term = dynamic_cast<RooAbsReal *>(iter->Next());
165 assert(uPdf != NULL); assert(extended_term != NULL); assert(iter->Next() == NULL);
166 delete iter;
167
168 RooAbsPdf *newUPdf = StripConstraints(*uPdf, observables);
169 if(newUPdf == NULL) return NULL; // only constraints in underlying pdf
170 else return new RooExtendPdf(TString::Format("%s_unconstrained", pdf.GetName()).Data(),
171 TString::Format("%s without constraints", pdf.GetTitle()).Data(), *newUPdf, *extended_term);
172
173 } else if (id == typeid(RooSimultaneous)) { //|| id == typeid(RooSimultaneousOpt)) {
174
175 RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf); assert(sim != NULL);
176 RooAbsCategoryLValue *cat = (RooAbsCategoryLValue *) sim->indexCat().Clone(); assert(cat != NULL);
177 RooArgList pdfList;
178
179 for (int ic = 0, nc = cat->numBins((const char *)NULL); ic < nc; ++ic) {
180 cat->setBin(ic);
181 RooAbsPdf* catPdf = sim->getPdf(cat->getCurrentLabel());
182 RooAbsPdf* newPdf = NULL;
183 // it is possible that a pdf is not defined for every category
184 if (catPdf != NULL) newPdf = StripConstraints(*catPdf, observables);
185 if (newPdf == NULL) { delete cat; return NULL; } // all channels must have observables
186 pdfList.add(*newPdf);
187 }
188
189 return new RooSimultaneous(TString::Format("%s_unconstrained", sim->GetName()).Data(),
190 TString::Format("%s without constraints", sim->GetTitle()).Data(), pdfList, *cat);
191
192 } else if (pdf.dependsOn(observables)) {
193 return (RooAbsPdf *) pdf.clone(TString::Format("%s_unconstrained", pdf.GetName()).Data());
194 }
195
196 return NULL; // just a constraint term
197 }
198
199 RooAbsPdf * MakeUnconstrainedPdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name) {
200 // make a clone pdf without all constraint terms in a common pdf
201 RooAbsPdf * unconstrainedPdf = StripConstraints(pdf, observables);
202 if(!unconstrainedPdf) {
203 oocoutE((TObject *)NULL, InputArguments) << "RooStats::MakeUnconstrainedPdf - invalid observable list passed (observables not found in original pdf) or invalid pdf passed (without observables)" << endl;
204 return NULL;
205 }
206 if(name != NULL) unconstrainedPdf->SetName(name);
207 return unconstrainedPdf;
208 }
209
211 // make a clone pdf without all constraint terms in a common pdf
212 if(!model.GetPdf() || !model.GetObservables()) {
213 oocoutE((TObject *)NULL, InputArguments) << "RooStatsUtils::MakeUnconstrainedPdf - invalid input model: missing pdf and/or observables" << endl;
214 return NULL;
215 }
216 return MakeUnconstrainedPdf(*model.GetPdf(), *model.GetObservables(), name);
217 }
218
219 // Helper class for GetAsTTree
220 class BranchStore {
221 public:
222 std::map<TString, Double_t> fVarVals;
223 double fInval;
224 TTree *fTree;
225
226 BranchStore(const vector <TString> &params = vector <TString>(), double _inval = -999.) : fTree(0) {
227 fInval = _inval;
228 for(unsigned int i = 0;i<params.size();i++)
229 fVarVals[params[i]] = _inval;
230 }
231
232 ~BranchStore() {
233 if (fTree) {
234 for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
235 TBranch *br = fTree->GetBranch( it->first );
236 if (br) br->ResetAddress();
237 }
238 }
239 }
240
241 void AssignToTTree(TTree &myTree) {
242 fTree = &myTree;
243 for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
244 const TString& name = it->first;
245 myTree.Branch( name, &fVarVals[name], TString::Format("%s/D", name.Data()));
246 }
247 }
248 void ResetValues() {
249 for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
250 const TString& name = it->first;
251 fVarVals[name] = fInval;
252 }
253 }
254 };
255
256 BranchStore* CreateBranchStore(const RooDataSet& data) {
257 if (data.numEntries() == 0) {
258 return new BranchStore;
259 }
260 vector <TString> V;
261 const RooArgSet* aset = data.get(0);
262 RooAbsArg *arg(0);
263 TIterator *it = aset->createIterator();
264 for(;(arg = dynamic_cast<RooAbsArg*>(it->Next()));) {
265 RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
266 if (rvar == NULL)
267 continue;
268 V.push_back(rvar->GetName());
269 if (rvar->hasAsymError()) {
270 V.push_back(TString::Format("%s_errlo", rvar->GetName()));
271 V.push_back(TString::Format("%s_errhi", rvar->GetName()));
272 }
273 else if (rvar->hasError()) {
274 V.push_back(TString::Format("%s_err", rvar->GetName()));
275 }
276 }
277 delete it;
278 return new BranchStore(V);
279 }
280
281 void FillTree(TTree &myTree, const RooDataSet &data) {
282 BranchStore *bs = CreateBranchStore(data);
283 bs->AssignToTTree(myTree);
284
285 for(int entry = 0;entry<data.numEntries();entry++) {
286 bs->ResetValues();
287 const RooArgSet* aset = data.get(entry);
288 RooAbsArg *arg(0);
289 RooLinkedListIter it = aset->iterator();
290 for(;(arg = dynamic_cast<RooAbsArg*>(it.Next()));) {
291 RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
292 if (rvar == NULL)
293 continue;
294 bs->fVarVals[rvar->GetName()] = rvar->getValV();
295 if (rvar->hasAsymError()) {
296 bs->fVarVals[TString::Format("%s_errlo", rvar->GetName())] = rvar->getAsymErrorLo();
297 bs->fVarVals[TString::Format("%s_errhi", rvar->GetName())] = rvar->getAsymErrorHi();
298 }
299 else if (rvar->hasError()) {
300 bs->fVarVals[TString::Format("%s_err", rvar->GetName())] = rvar->getError();
301 }
302 }
303 myTree.Fill();
304 }
305
306 delete bs;
307 }
308
310 TTree* myTree = new TTree(name, desc);
311 FillTree(*myTree, data);
312 return myTree;
313 }
314
315
316 // useful function to print in one line the content of a set with their values
317 void PrintListContent(const RooArgList & l, std::ostream & os ) {
318 bool first = true;
319 os << "( ";
320 for (int i = 0; i< l.getSize(); ++i) {
321 if (first) {
322 first=kFALSE ;
323 } else {
324 os << ", " ;
325 }
326 l[i].printName(os);
327 os << " = ";
328 l[i].printValue(os);
329 }
330 os << ")\n";
331 }
332}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define oocoutW(o, a)
Definition: RooMsgService.h:47
#define oocoutE(o, a)
Definition: RooMsgService.h:48
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
#define NamespaceImp(name)
Definition: Rtypes.h:376
char name[80]
Definition: TGX11.cxx:109
double sqrt(double)
double log(double)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
Bool_t dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, Bool_t valueOnly=kFALSE) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:730
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:85
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
virtual TObject * clone(const char *newname=0) const =0
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2216
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
virtual void setBin(Int_t ibin, const char *rangeName=0)
Set category to i-th fit bin, which is the i-th registered state.
virtual Int_t numBins(const char *rangeName) const
Return the number of fit bins ( = number of types )
virtual const char * getCurrentLabel() const
Return label string of current state.
Int_t getSize() const
Bool_t contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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'.
RooExtendPdf is a wrapper around an existing PDF that adds a parameteric extended likelihood term to ...
Definition: RooExtendPdf.h:22
A wrapper around TIterator derivatives.
TObject * Next() override
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
const RooArgList & pdfList() const
Definition: RooProdPdf.h:67
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
Double_t getAsymErrorLo() const
Definition: RooRealVar.h:66
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:68
Double_t getAsymErrorHi() const
Definition: RooRealVar.h:67
virtual Double_t getValV(const RooArgSet *nset=0) const
Return value of variable.
Definition: RooRealVar.cxx:250
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:63
Double_t getError() const
Definition: RooRealVar.h:62
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
const RooAbsCategoryLValue & indexCat() const
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:249
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
A TTree is a list of TBranches.
Definition: TBranch.h:91
virtual void ResetAddress()
Reset the address of the branch.
Definition: TBranch.cxx:2515
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
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
Mother of all ROOT objects.
Definition: TObject.h:37
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2311
A TTree represents a columnar dataset.
Definition: TTree.h:78
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4524
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:5209
TBranch * Branch(const char *name, T *obj, Int_t bufsize=32000, Int_t splitlevel=99)
Add a new branch, and infer the data type from the type of obj being passed.
Definition: TTree.h:348
const Int_t n
Definition: legend1.C:16
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:19
Double_t AsimovSignificance(Double_t s, Double_t b, Double_t sigma_b=0.0)
Compute the Asimov Median significance for a Poisson process with s = expected number of signal event...
TTree * GetAsTTree(TString name, TString desc, const RooDataSet &data)
RooAbsPdf * MakeUnconstrainedPdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name=NULL)
void FillTree(TTree &myTree, const RooDataSet &data)
BranchStore * CreateBranchStore(const RooDataSet &data)
RooAbsPdf * StripConstraints(RooAbsPdf &pdf, const RooArgSet &observables)
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
void FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints)
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
void UseNLLOffset(bool on)
Use an offset in NLL calculations.
bool IsNLLOffset()
Test of RooStats should by default offset NLL calculations.
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
static constexpr double s
Definition: first.py:1
auto * l
Definition: textangle.C:4