Logo ROOT   master
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"
25 #include "RooStats/ModelConfig.h"
26 #include "RooStats/RooStatsUtils.h"
27 #include <typeinfo>
28 
29 using namespace std;
30 
31 namespace 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  } else if(constraints.getSize() == 1) {
126  return dynamic_cast<RooAbsPdf *>(constraints.first()->clone(name));
127  }
128  return new RooProdPdf(name,"", constraints);
129  }
130 
131  RooAbsPdf * MakeNuisancePdf(const RooStats::ModelConfig &model, const char *name) {
132  // make a nuisance pdf by factorizing out all constraint terms in a common pdf
133  if (!model.GetPdf() || !model.GetObservables() ) {
134  oocoutE((TObject*)0, InputArguments) << "RooStatsUtils::MakeNuisancePdf - invalid input model: missing pdf and/or observables" << endl;
135  return 0;
136  }
137  return MakeNuisancePdf(*model.GetPdf(), *model.GetObservables(), name);
138  }
139 
140  RooAbsPdf * StripConstraints(RooAbsPdf &pdf, const RooArgSet &observables) {
141  const std::type_info & id = typeid(pdf);
142 
143  if (id == typeid(RooProdPdf)) {
144 
145  RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
146  RooArgList list(prod->pdfList()); RooArgList newList;
147 
148  for (int i = 0, n = list.getSize(); i < n; ++i) {
149  RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
150  RooAbsPdf *newPdfi = StripConstraints(*pdfi, observables);
151  if(newPdfi != NULL) newList.add(*newPdfi);
152  }
153 
154  if(newList.getSize() == 0) return NULL; // only constraints in product
155  // return single component (no longer a product)
156  else if(newList.getSize() == 1) return dynamic_cast<RooAbsPdf *>(newList.at(0)->clone(TString::Format("%s_unconstrained",
157  newList.at(0)->GetName())));
158  else return new RooProdPdf(TString::Format("%s_unconstrained", prod->GetName()).Data(),
159  TString::Format("%s without constraints", prod->GetTitle()).Data(), newList);
160 
161  } else if (id == typeid(RooExtendPdf)) {
162 
163  TIterator *iter = pdf.serverIterator();
164  // extract underlying pdf which is extended; first server is the pdf; second server is the number of events variable
165  RooAbsPdf *uPdf = dynamic_cast<RooAbsPdf *>(iter->Next());
166  RooAbsReal *extended_term = dynamic_cast<RooAbsReal *>(iter->Next());
167  assert(uPdf != NULL); assert(extended_term != NULL); assert(iter->Next() == NULL);
168  delete iter;
169 
170  RooAbsPdf *newUPdf = StripConstraints(*uPdf, observables);
171  if(newUPdf == NULL) return NULL; // only constraints in underlying pdf
172  else return new RooExtendPdf(TString::Format("%s_unconstrained", pdf.GetName()).Data(),
173  TString::Format("%s without constraints", pdf.GetTitle()).Data(), *newUPdf, *extended_term);
174 
175  } else if (id == typeid(RooSimultaneous)) { //|| id == typeid(RooSimultaneousOpt)) {
176 
177  RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf); assert(sim != NULL);
178  RooAbsCategoryLValue *cat = (RooAbsCategoryLValue *) sim->indexCat().Clone(); assert(cat != NULL);
179  RooArgList pdfList;
180 
181  for (int ic = 0, nc = cat->numBins((const char *)NULL); ic < nc; ++ic) {
182  cat->setBin(ic);
183  RooAbsPdf* catPdf = sim->getPdf(cat->getCurrentLabel());
184  RooAbsPdf* newPdf = NULL;
185  // it is possible that a pdf is not defined for every category
186  if (catPdf != NULL) newPdf = StripConstraints(*catPdf, observables);
187  if (newPdf == NULL) { delete cat; return NULL; } // all channels must have observables
188  pdfList.add(*newPdf);
189  }
190 
191  return new RooSimultaneous(TString::Format("%s_unconstrained", sim->GetName()).Data(),
192  TString::Format("%s without constraints", sim->GetTitle()).Data(), pdfList, *cat);
193 
194  } else if (pdf.dependsOn(observables)) {
195  return (RooAbsPdf *) pdf.clone(TString::Format("%s_unconstrained", pdf.GetName()).Data());
196  }
197 
198  return NULL; // just a constraint term
199  }
200 
201  RooAbsPdf * MakeUnconstrainedPdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name) {
202  // make a clone pdf without all constraint terms in a common pdf
203  RooAbsPdf * unconstrainedPdf = StripConstraints(pdf, observables);
204  if(!unconstrainedPdf) {
205  oocoutE((TObject *)NULL, InputArguments) << "RooStats::MakeUnconstrainedPdf - invalid observable list passed (observables not found in original pdf) or invalid pdf passed (without observables)" << endl;
206  return NULL;
207  }
208  if(name != NULL) unconstrainedPdf->SetName(name);
209  return unconstrainedPdf;
210  }
211 
213  // make a clone pdf without all constraint terms in a common pdf
214  if(!model.GetPdf() || !model.GetObservables()) {
215  oocoutE((TObject *)NULL, InputArguments) << "RooStatsUtils::MakeUnconstrainedPdf - invalid input model: missing pdf and/or observables" << endl;
216  return NULL;
217  }
218  return MakeUnconstrainedPdf(*model.GetPdf(), *model.GetObservables(), name);
219  }
220 
221  // Helper class for GetAsTTree
222  class BranchStore {
223  public:
224  std::map<TString, Double_t> fVarVals;
225  double fInval;
226  TTree *fTree;
227 
228  BranchStore(const vector <TString> &params = vector <TString>(), double _inval = -999.) : fTree(0) {
229  fInval = _inval;
230  for(unsigned int i = 0;i<params.size();i++)
231  fVarVals[params[i]] = _inval;
232  }
233 
234  ~BranchStore() {
235  if (fTree) {
236  for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
237  TBranch *br = fTree->GetBranch( it->first );
238  if (br) br->ResetAddress();
239  }
240  }
241  }
242 
243  void AssignToTTree(TTree &myTree) {
244  fTree = &myTree;
245  for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
246  const TString& name = it->first;
247  myTree.Branch( name, &fVarVals[name], TString::Format("%s/D", name.Data()));
248  }
249  }
250  void ResetValues() {
251  for(std::map<TString, Double_t>::iterator it = fVarVals.begin();it!=fVarVals.end();++it) {
252  const TString& name = it->first;
253  fVarVals[name] = fInval;
254  }
255  }
256  };
257 
258  BranchStore* CreateBranchStore(const RooDataSet& data) {
259  if (data.numEntries() == 0) {
260  return new BranchStore;
261  }
262  vector <TString> V;
263  const RooArgSet* aset = data.get(0);
264  RooAbsArg *arg(0);
265  TIterator *it = aset->createIterator();
266  for(;(arg = dynamic_cast<RooAbsArg*>(it->Next()));) {
267  RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
268  if (rvar == NULL)
269  continue;
270  V.push_back(rvar->GetName());
271  if (rvar->hasAsymError()) {
272  V.push_back(TString::Format("%s_errlo", rvar->GetName()));
273  V.push_back(TString::Format("%s_errhi", rvar->GetName()));
274  }
275  else if (rvar->hasError()) {
276  V.push_back(TString::Format("%s_err", rvar->GetName()));
277  }
278  }
279  delete it;
280  return new BranchStore(V);
281  }
282 
283  void FillTree(TTree &myTree, const RooDataSet &data) {
284  BranchStore *bs = CreateBranchStore(data);
285  bs->AssignToTTree(myTree);
286 
287  for(int entry = 0;entry<data.numEntries();entry++) {
288  bs->ResetValues();
289  const RooArgSet* aset = data.get(entry);
290  RooAbsArg *arg(0);
291  RooLinkedListIter it = aset->iterator();
292  for(;(arg = dynamic_cast<RooAbsArg*>(it.Next()));) {
293  RooRealVar *rvar = dynamic_cast<RooRealVar*>(arg);
294  if (rvar == NULL)
295  continue;
296  bs->fVarVals[rvar->GetName()] = rvar->getValV();
297  if (rvar->hasAsymError()) {
298  bs->fVarVals[TString::Format("%s_errlo", rvar->GetName())] = rvar->getAsymErrorLo();
299  bs->fVarVals[TString::Format("%s_errhi", rvar->GetName())] = rvar->getAsymErrorHi();
300  }
301  else if (rvar->hasError()) {
302  bs->fVarVals[TString::Format("%s_err", rvar->GetName())] = rvar->getError();
303  }
304  }
305  myTree.Fill();
306  }
307 
308  delete bs;
309  }
310 
311  TTree * GetAsTTree(TString name, TString desc, const RooDataSet& data) {
312  TTree* myTree = new TTree(name, desc);
313  FillTree(*myTree, data);
314  return myTree;
315  }
316 
317 
318  // useful function to print in one line the content of a set with their values
319  void PrintListContent(const RooArgList & l, std::ostream & os ) {
320  bool first = true;
321  os << "( ";
322  for (int i = 0; i< l.getSize(); ++i) {
323  if (first) {
324  first=kFALSE ;
325  } else {
326  os << ", " ;
327  }
328  l[i].printName(os);
329  os << " = ";
330  l[i].printValue(os);
331  }
332  os << ")\n";
333  }
334 }
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual TObject * clone(const char *newname=0) const =0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
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
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:731
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooProdPdf is an efficient implementation of a product of PDFs of the form .
Definition: RooProdPdf.h:31
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4507
Basic string class.
Definition: TString.h:131
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:85
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...
STL namespace.
virtual Double_t getValV(const RooArgSet *nset=0) const
Return value of variable.
Definition: RooRealVar.cxx:214
TTree * GetAsTTree(TString name, TString desc, const RooDataSet &data)
Iterator abstract base class.
Definition: TIterator.h:30
TObject * Next() override
void FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints)
double sqrt(double)
virtual Int_t numBins(const char *rangeName) const
Return the number of fit bins ( = number of types )
RooAbsPdf * StripConstraints(RooAbsPdf &pdf, const RooArgSet &observables)
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:68
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
#define oocoutE(o, a)
Definition: RooMsgService.h:48
virtual const char * getCurrentLabel() const
Return label string of current state.
static constexpr double s
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:63
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:5163
Int_t getSize() const
ROOT::R::TRInterface & r
Definition: Object.C:4
RooAbsPdf * MakeUnconstrainedPdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name=NULL)
Double_t getAsymErrorHi() const
Definition: RooRealVar.h:67
const RooAbsCategoryLValue & indexCat() const
RooAbsArg * first() const
const RooArgList & pdfList() const
Definition: RooProdPdf.h:67
BranchStore * CreateBranchStore(const RooDataSet &data)
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2217
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
#define NamespaceImp(name)
Definition: Rtypes.h:376
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
void UseNLLOffset(bool on)
Use an offset in NLL calculations.
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
const Bool_t kFALSE
Definition: RtypesCore.h:90
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:345
Namespace for the RooStats classes.
Definition: Asimov.h:19
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
virtual void ResetAddress()
Reset the address of the branch.
Definition: TBranch.cxx:2465
double Double_t
Definition: RtypesCore.h:57
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
bool IsNLLOffset()
Test of RooStats should by default offset NLL calculations.
#define oocoutW(o, a)
Definition: RooMsgService.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
RooExtendPdf is a wrapper around an existing PDF that adds a parameteric extended likelihood term to ...
Definition: RooExtendPdf.h:22
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
void FillTree(TTree &myTree, const RooDataSet &data)
auto * l
Definition: textangle.C:4
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
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 TObject * Next()=0
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
A TTree represents a columnar dataset.
Definition: TTree.h:78
Bool_t contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
Definition: first.py:1
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
A TTree is a list of TBranches.
Definition: TBranch.h:88
Double_t getError() const
Definition: RooRealVar.h:62
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Int_t n
Definition: legend1.C:16
Double_t getAsymErrorLo() const
Definition: RooRealVar.h:66
A wrapper around TIterator derivatives.
char name[80]
Definition: TGX11.cxx:109
double log(double)
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48