ROOT   Reference Guide
Searching...
No Matches
RooStatsUtils.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer 28/07/2008
3
4/*************************************************************************
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)) {
101 } else {
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);
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
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 );
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);
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)
#define oocoutE(o, a)
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
#define NamespaceImp(name)
Definition Rtypes.h:379
char name[80]
Definition TGX11.cxx:110
double sqrt(double)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:72
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.
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition RooAbsArg.h:84
TIterator * serverIterator() const
Definition RooAbsArg.h:140
virtual TObject * clone(const char *newname=0) const =0
void SetName(const char *name)
Set the name of the TNamed.
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.
TIterator-style iteration over contained elements.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:61
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:70
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
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 ...
A wrapper around TIterator derivatives.
TObject * Next() override
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:37
const RooArgList & pdfList() const
Definition RooProdPdf.h:73
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
Double_t getAsymErrorLo() const
Definition RooRealVar.h:64
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition RooRealVar.h:66
Double_t getAsymErrorHi() const
Definition RooRealVar.h:65
virtual Double_t getValV(const RooArgSet *nset=0) const
Return value of variable.
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition RooRealVar.h:61
Double_t getError() const
Definition RooRealVar.h:60
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)
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
A TTree is a list of TBranches.
Definition TBranch.h:89
Reset the address of the branch.
Definition TBranch.cxx:2589
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:136
const char * Data() const
Definition TString.h:369
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:2331
A TTree represents a columnar dataset.
Definition TTree.h:79
virtual Int_t Fill()
Fill all branches.
Definition TTree.cxx:4590
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition TTree.cxx:5275
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:350
const Int_t n
Definition legend1.C:16
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)
Definition first.py:1
auto * l
Definition textangle.C:4