Logo ROOT   6.14/05
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"
24 #include "RooStats/ModelConfig.h"
25 #include "RooStats/RooStatsUtils.h"
26 #include <typeinfo>
27 
28 using namespace std;
29 
30 // this file is only for the documentation of RooStats namespace
31 
32 namespace 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 
105  void FactorizePdf(RooStats::ModelConfig &model, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
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 }
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
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:237
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:735
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
virtual void printValue(std::ostream &os) const
Print value of collection, i.e.
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
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:193
TTree * GetAsTTree(TString name, TString desc, const RooDataSet &data)
Iterator abstract base class.
Definition: TIterator.h:30
void FactorizePdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints)
double sqrt(double)
virtual Int_t numBins(const char *rangeName) const
Returm 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:59
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
#define oocoutE(o, a)
Definition: RooMsgService.h:47
virtual TObject * Next()
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:54
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void printName(std::ostream &os) const
Return collection name.
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:58
const RooAbsCategoryLValue & indexCat() const
RooAbsArg * first() const
const RooArgList & pdfList() const
Definition: RooProdPdf.h:69
BranchStore * CreateBranchStore(const RooDataSet &data)
void SetName(const char *name)
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2381
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
#define NamespaceImp(name)
Definition: Rtypes.h:374
virtual const char * getLabel() const
Return label string of current state.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
void UseNLLOffset(bool on)
TIterator * serverIterator() const
Definition: RooAbsArg.h:112
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
const Bool_t kFALSE
Definition: RtypesCore.h:88
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:222
virtual void ResetAddress()
Reset the address of the branch.
Definition: TBranch.cxx:2227
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
bool gUseOffset
static constexpr double s
bool IsNLLOffset()
#define oocoutW(o, a)
Definition: RooMsgService.h:46
Mother of all ROOT objects.
Definition: TObject.h:37
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
void FillTree(TTree &myTree, const RooDataSet &data)
auto * l
Definition: textangle.C:4
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
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
Bool_t contains(const RooAbsArg &var) const
Definition: first.py:1
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
A TTree is a list of TBranches.
Definition: TBranch.h:62
Double_t getError() const
Definition: RooRealVar.h:53
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:57
RooLinkedListIter is the TIterator implementation for RooLinkedList.
char name[80]
Definition: TGX11.cxx:109
double log(double)
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:285
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48