Logo ROOT   6.16/01
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
20#include "RooUniform.h"
21#include "RooProdPdf.h"
22#include "RooExtendPdf.h"
23#include "RooSimultaneous.h"
26#include <typeinfo>
27
28using namespace std;
29
30// this file is only for the documentation of RooStats namespace
31
32namespace RooStats {
33
34 bool gUseOffset = false;
35
37 // Asimov significance
38 // formula [10] and [20] from https://www.pp.rhul.ac.uk/~cowan/stat/notes/medsigNote.pdf
39 // case we have a sigma_b
40 double sb2 = sigma_b*sigma_b;
41 // formula below has a large error when sigma_b becomes zero
42 // better to use the approximation for sigma_b=0 for very small values
43 double r = sb2/b;
44 if (r > 1.E-12) {
45 double bpsb2 = b + sb2;
46 double b2 = b*b;
47 double spb = s+b;
48 double za2 = 2.*( (spb)* std::log( ( spb)*(bpsb2)/(b2+ spb*sb2) ) -
49 (b2/sb2) * std::log(1. + ( sb2 * s)/(b * bpsb2) ) );
50 return sqrt(za2);
51
52 }
53 // case when the background (b) is known
54 double za2 = 2.*( (s+b) * std::log(1. + s/b) -s );
55 return std::sqrt(za2);
56 }
57
58
59 void UseNLLOffset(bool on) {
60 // use offset in NLL calculations
61 gUseOffset = on;
62 }
63
64 bool IsNLLOffset() {
65 return gUseOffset;
66 }
67
68 void FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
69 // utility function to factorize constraint terms from a pdf
70 // (from G. Petrucciani)
71 const std::type_info & id = typeid(pdf);
72 if (id == typeid(RooProdPdf)) {
73 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
74 RooArgList list(prod->pdfList());
75 for (int i = 0, n = list.getSize(); i < n; ++i) {
76 RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
77 FactorizePdf(observables, *pdfi, obsTerms, constraints);
78 }
79 } else if (id == typeid(RooExtendPdf)) {
80 TIterator *iter = pdf.serverIterator();
81 // extract underlying pdf which is extended; first server is the pdf; second server is the number of events variable
82 RooAbsPdf *updf = dynamic_cast<RooAbsPdf *>(iter->Next());
83 assert(updf != 0);
84 delete iter;
85 FactorizePdf(observables, *updf, obsTerms, constraints);
86 } else if (id == typeid(RooSimultaneous)) { //|| id == typeid(RooSimultaneousOpt)) {
87 RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf);
88 assert(sim != 0);
90 for (int ic = 0, nc = cat->numBins((const char *)0); ic < nc; ++ic) {
91 cat->setBin(ic);
92 RooAbsPdf* catPdf = sim->getPdf(cat->getLabel());
93 // it is possible that a pdf is not defined for every category
94 if (catPdf != 0) FactorizePdf(observables, *catPdf, obsTerms, constraints);
95 }
96 delete cat;
97 } else if (pdf.dependsOn(observables)) {
98 if (!obsTerms.contains(pdf)) obsTerms.add(pdf);
99 } else {
100 if (!constraints.contains(pdf)) constraints.add(pdf);
101 }
102 }
103
104
106 // utility function to factorize constraint terms from a pdf
107 // (from G. Petrucciani)
108 if (!model.GetObservables() ) {
109 oocoutE((TObject*)0,InputArguments) << "RooStatsUtils::FactorizePdf - invalid input model: missing observables" << endl;
110 return;
111 }
112 return FactorizePdf(*model.GetObservables(), pdf, obsTerms, constraints);
113 }
114
115
116 RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name) {
117 // make a nuisance pdf by factorizing out all constraint terms in a common pdf
118 RooArgList obsTerms, constraints;
119 FactorizePdf(observables, pdf, obsTerms, constraints);
120 if(constraints.getSize() == 0) {
121 oocoutW((TObject *)0, Eval) << "RooStatsUtils::MakeNuisancePdf - no constraints found on nuisance parameters in the input model" << endl;
122 return 0;
123 } else if(constraints.getSize() == 1) {
124 return dynamic_cast<RooAbsPdf *>(constraints.first()->clone(name));
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->getLabel());
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
309 TTree * GetAsTTree(TString name, TString desc, const RooDataSet& data) {
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:46
#define oocoutE(o, a)
Definition: RooMsgService.h:47
const Bool_t kFALSE
Definition: RtypesCore.h:88
double Double_t
Definition: RtypesCore.h:55
#define NamespaceImp(name)
Definition: Rtypes.h:378
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:66
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:736
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
TIterator * serverIterator() const
Definition: RooAbsArg.h:112
virtual TObject * clone(const char *newname=0) const =0
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2382
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
Returm the number of fit bins ( = number of types )
virtual const char * getLabel() const
Return label string of current state.
Int_t getSize() const
Bool_t contains(const RooAbsArg &var) const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
RooAbsArg * first() const
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator * createIterator(Bool_t dir=kIterForward) const
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
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:31
RooExtendPdf is a wrapper around an existing PDF that adds a parameteric extended likelihood term to ...
Definition: RooExtendPdf.h:22
RooLinkedListIter is the TIterator implementation for RooLinkedList.
virtual TObject * Next()
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
const RooArgList & pdfList() const
Definition: RooProdPdf.h:69
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
Double_t getAsymErrorLo() const
Definition: RooRealVar.h:57
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:59
Double_t getAsymErrorHi() const
Definition: RooRealVar.h:58
virtual Double_t getValV(const RooArgSet *nset=0) const
Return value of variable.
Definition: RooRealVar.cxx:193
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:54
Double_t getError() const
Definition: RooRealVar.h:53
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
A TTree is a list of TBranches.
Definition: TBranch.h:64
virtual void ResetAddress()
Reset the address of the branch.
Definition: TBranch.cxx:2234
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
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:2286
const Int_t n
Definition: legend1.C:16
@ InputArguments
Definition: RooGlobalFunc.h:58
@(#)root/roostats:$Id$ Author: George Lewis, Kyle Cranmer
Definition: Asimov.h:20
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)
void UseNLLOffset(bool on)
bool IsNLLOffset()
bool gUseOffset
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
static constexpr double s
Definition: first.py:1
STL namespace.
auto * l
Definition: textangle.C:4