Logo ROOT   6.10/09
Reference Guide
HistFactoryModelUtils.cxx
Go to the documentation of this file.
1 
2 // A set of utils for navegating HistFactory models
3 #include <stdexcept>
4 #include <typeinfo>
5 
7 #include "TIterator.h"
8 #include "RooAbsArg.h"
9 #include "RooAbsPdf.h"
10 #include "RooArgSet.h"
11 #include "RooArgList.h"
12 #include "RooSimultaneous.h"
13 #include "RooCategory.h"
14 #include "RooRealVar.h"
15 #include "RooProdPdf.h"
16 #include "TH1.h"
17 
20 
21 namespace RooStats{
22 namespace HistFactory{
23 
24 
25  std::string channelNameFromPdf( RooAbsPdf* channelPdf ) {
26  std::string channelPdfName = channelPdf->GetName();
27  std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
28  return ChannelName;
29  }
30 
32 
33  bool verbose=false;
34 
35  if(verbose) std::cout << "Getting the RooRealSumPdf for the channel: "
36  << sim_channel->GetName() << std::endl;
37 
38  std::string channelPdfName = sim_channel->GetName();
39  std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
40 
41  // Now, get the RooRealSumPdf
42  // ie the channel WITHOUT constraints
43  std::string realSumPdfName = ChannelName + "_model";
44 
45  RooAbsPdf* sum_pdf = NULL;
46  TIterator* iter_sum_pdf = sim_channel->getComponents()->createIterator(); //serverIterator();
47  bool FoundSumPdf=false;
48  RooAbsArg* sum_pdf_arg=NULL;
49  while((sum_pdf_arg=(RooAbsArg*)iter_sum_pdf->Next())) {
50  std::string NodeClassName = sum_pdf_arg->ClassName();
51  if( NodeClassName == std::string("RooRealSumPdf") ) {
52  FoundSumPdf=true;
53  sum_pdf = (RooAbsPdf*) sum_pdf_arg;
54  break;
55  }
56  }
57  if( ! FoundSumPdf ) {
58  if(verbose) {
59  std::cout << "Failed to find RooRealSumPdf for channel: " << sim_channel->GetName() << std::endl;
60  sim_channel->getComponents()->Print("V");
61  }
62  sum_pdf=NULL;
63  //throw std::runtime_error("Failed to find RooRealSumPdf for channel");
64  }
65  else {
66  if(verbose) std::cout << "Found RooRealSumPdf: " << sum_pdf->GetName() << std::endl;
67  }
68  delete iter_sum_pdf;
69  iter_sum_pdf = NULL;
70 
71  return sum_pdf;
72 
73  }
74 
75 
76  void FactorizeHistFactoryPdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
77  // utility function to factorize constraint terms from a pdf
78  // (from G. Petrucciani)
79  const std::type_info & id = typeid(pdf);
80  if (id == typeid(RooProdPdf)) {
81  RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
82  RooArgList list(prod->pdfList());
83  for (int i = 0, n = list.getSize(); i < n; ++i) {
84  RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
85  FactorizeHistFactoryPdf(observables, *pdfi, obsTerms, constraints);
86  }
87  } else if (id == typeid(RooSimultaneous) || id == typeid(HistFactorySimultaneous) ) { //|| id == typeid(RooSimultaneousOpt)) {
88  RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf);
90  for (int ic = 0, nc = cat->numBins((const char *)0); ic < nc; ++ic) {
91  cat->setBin(ic);
92  FactorizeHistFactoryPdf(observables, *sim->getPdf(cat->getLabel()), obsTerms, constraints);
93  }
94  delete cat;
95  } else if (pdf.dependsOn(observables)) {
96  if (!obsTerms.contains(pdf)) obsTerms.add(pdf);
97  } else {
98  if (!constraints.contains(pdf)) constraints.add(pdf);
99  }
100  }
101 
102  /*
103  void getChannelsFromModel( RooAbsPdf* model, RooArgSet* channels, RooArgSet* channelsWithConstraints ) {
104 
105  // Loop through the model
106  // Find all channels
107 
108  std::string modelClassName = model->ClassName();
109 
110  if( modelClassName == std::string("RooSimultaneous") || model->InheritsFrom("RooSimultaneous") ) {
111 
112  TIterator* simServerItr = model->serverIterator();
113 
114  // Loop through the child nodes of the sim pdf
115  // and find the channel nodes
116  RooAbsArg* sim_channel_arg = NULL;
117  while(( sim_channel = (RooAbsArg*) simServerItr->Next() )) {
118 
119  RooAbsPdf* sim_channel = (RooAbsPdf*) sim_channel_arg;
120 
121  // Ignore the Channel Cat
122  std::string channelPdfName = sim_channel->GetName();
123  std::string channelClassName = sim_channel->ClassName();
124  if( channelClassName == std::string("RooCategory") ) continue;
125 
126  // If we got here, we found a channel.
127  // Format is model_<ChannelName>
128 
129  std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
130 
131  // Now, get the RooRealSumPdf
132  RooAbsPdf* sum_pdf = getSumPdfFromChannel( sim_channel );
133 
134 
135  / *
136  // Now, get the RooRealSumPdf
137  // ie the channel WITHOUT constraints
138 
139  std::string realSumPdfName = ChannelName + "_model";
140 
141  RooAbsPdf* sum_pdf = NULL;
142  TIterator* iter_sum_pdf = sim_channel->getComponents()->createIterator(); //serverIterator();
143  bool FoundSumPdf=false;
144  RooAbsArg* sum_pdf_arg=NULL;
145  while((sum_pdf_arg=(RooAbsArg*)iter_sum_pdf->Next())) {
146 
147  std::string NodeClassName = sum_pdf_arg->ClassName();
148  if( NodeClassName == std::string("RooRealSumPdf") ) {
149  FoundSumPdf=true;
150  sum_pdf = (RooAbsPdf*) sum_pdf_arg;
151  break;
152  }
153  }
154  if( ! FoundSumPdf ) {
155  std::cout << "Failed to find RooRealSumPdf for channel: " << sim_channel->GetName() << std::endl;
156  sim_channel->getComponents()->Print("V");
157  throw std::runtime_error("Failed to find RooRealSumPdf for channel");
158  }
159  delete iter_sum_pdf;
160  iter_sum_pdf = NULL;
161  * /
162 
163  // Okay, now add to the arg sets
164  channels->add( *sum_pdf );
165  channelsWithConstraints->add( *sim_channel );
166 
167  }
168 
169  delete simServerItr;
170 
171  }
172  else {
173  std::cout << "Model is not a RooSimultaneous or doesn't derive from one." << std::endl;
174  std::cout << "HistFactoryModelUtils isn't yet implemented for these pdf's" << std::endl;
175  }
176 
177  }
178  */
179 
180  bool getStatUncertaintyFromChannel( RooAbsPdf* channel, ParamHistFunc*& paramfunc, RooArgList* gammaList ) {
181 
182  bool verbose=false;
183 
184  // Find the servers of this channel
185  //TIterator* iter = channel->serverIterator();
186  TIterator* iter = channel->getComponents()->createIterator(); //serverIterator();
187  bool FoundParamHistFunc=false;
188  RooAbsArg* paramfunc_arg = NULL;
189  while(( paramfunc_arg = (RooAbsArg*) iter->Next() )) {
190  std::string NodeName = paramfunc_arg->GetName();
191  std::string NodeClassName = paramfunc_arg->ClassName();
192  if( NodeClassName != std::string("ParamHistFunc") ) continue;
193  if( NodeName.find("mc_stat_") != std::string::npos ) {
194  FoundParamHistFunc=true;
195  paramfunc = (ParamHistFunc*) paramfunc_arg;
196  break;
197  }
198  }
199  if( ! FoundParamHistFunc || !paramfunc ) {
200  if(verbose) std::cout << "Failed to find ParamHistFunc for channel: " << channel->GetName() << std::endl;
201  return false;
202  }
203 
204  delete iter;
205  iter = NULL;
206 
207  // Now, get the set of gamma's
208  gammaList = (RooArgList*) &( paramfunc->paramList());
209  if(verbose) gammaList->Print("V");
210 
211  return true;
212 
213  }
214 
215 
216  void getDataValuesForObservables( std::map< std::string, std::vector<double> >& ChannelBinDataMap,
217  RooAbsData* data, RooAbsPdf* pdf ) {
218 
219  bool verbose=false;
220 
221  //std::map< std::string, std::vector<int> ChannelBinDataMap;
222 
223  RooSimultaneous* simPdf = (RooSimultaneous*) pdf;
224 
225  // get category label
226  RooArgSet* allobs = (RooArgSet*) data->get();
227  TIterator* obsIter = allobs->createIterator();
228  RooCategory* cat = NULL;
229  RooAbsArg* temp = NULL;
230  while( (temp=(RooAbsArg*) obsIter->Next())) {
231  // use dynamic cast here instead
232  if( strcmp(temp->ClassName(),"RooCategory")==0){
233  cat = (RooCategory*) temp;
234  break;
235  }
236  }
237  if(verbose) {
238  if(!cat) std::cout <<"didn't find category"<< std::endl;
239  else std::cout <<"found category"<< std::endl;
240  }
241  delete obsIter;
242 
243  // split dataset
244  TList* dataByCategory = data->split(*cat);
245  if(verbose) dataByCategory->Print();
246  // note :
247  // RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject("");
248 
249  // loop over channels
250  RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
251  TIterator* iter = channelCat->typeIterator() ;
252  RooCatType* tt = NULL;
253  while((tt=(RooCatType*) iter->Next())) {
254 
255  // Get pdf associated with state from simpdf
256  RooAbsPdf* pdftmp = simPdf->getPdf(tt->GetName()) ;
257 
258  std::string ChannelName = pdftmp->GetName(); //tt->GetName();
259  if(verbose) std::cout << "Getting data for channel: " << ChannelName << std::endl;
260  ChannelBinDataMap[ ChannelName ] = std::vector<double>();
261 
262  RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject(tt->GetName());
263  if(verbose) dataForChan->Print();
264 
265  // Generate observables defined by the pdf associated with this state
266  RooArgSet* obstmp = pdftmp->getObservables(*dataForChan->get()) ;
267  RooRealVar* obs = ((RooRealVar*)obstmp->first());
268  if(verbose) obs->Print();
269 
270  //double expected = pdftmp->expectedEvents(*obstmp);
271 
272  // set value to desired value (this is just an example)
273  // double obsVal = obs->getVal();
274  // set obs to desired value of observable
275  // obs->setVal( obsVal );
276  //double fracAtObsValue = pdftmp->getVal(*obstmp);
277 
278  // get num events expected in bin for obsVal
279  // double nu = expected * fracAtObsValue;
280 
281  // an easier way to get n
282  TH1* histForN = dataForChan->createHistogram("HhstForN",*obs);
283  for(int i=1; i<=histForN->GetNbinsX(); ++i){
284  double n = histForN->GetBinContent(i);
285  if(verbose) std::cout << "n" << i << " = " << n << std::endl;
286  ChannelBinDataMap[ ChannelName ].push_back( n );
287  }
288  delete histForN;
289 
290  } // End Loop Over Categories
291 
292  delete iter;
293  return;
294 
295  }
296 
297 
298  int getStatUncertaintyConstraintTerm( RooArgList* constraints, RooRealVar* gamma_stat,
299  RooAbsReal*& pois_nom, RooRealVar*& tau ) {
300  // Given a set of constraint terms,
301  // find the poisson constraint for the
302  // given gamma and return the mean
303  // as well as the 'tau' parameter
304 
305  bool verbose=false;
306 
307  // To get the constraint term, loop over all constraint terms
308  // and look for the gamma_stat name as well as '_constraint'
309  // std::string constraintTermName = std::string(gamma_stat->GetName()) + "_constraint";
310  TIterator* iter_list = constraints->createIterator();
311  RooAbsArg* term_constr=NULL;
312  bool FoundConstraintTerm=false;
313  RooAbsPdf* constraintTerm=NULL;
314  while((term_constr=(RooAbsArg*)iter_list->Next())) {
315  std::string TermName = term_constr->GetName();
316  // std::cout << "Checking if is a constraint term: " << TermName << std::endl;
317 
318  //if( TermName.find(gamma_stat->GetName())!=string::npos ) {
319  if( term_constr->dependsOn( *gamma_stat) ) {
320  if( TermName.find("_constraint")!=std::string::npos ) {
321  FoundConstraintTerm=true;
322  constraintTerm = (RooAbsPdf*) term_constr;
323  break;
324  }
325  }
326  }
327  if( FoundConstraintTerm==false ) {
328  std::cout << "Error: Couldn't find constraint term for parameter: " << gamma_stat->GetName()
329  << " among constraints: " << constraints->GetName() << std::endl;
330  constraints->Print("V");
331  throw std::runtime_error("Failed to find Gamma ConstraintTerm");
332  return -1;
333  }
334  delete iter_list;
335 
336  /*
337  RooAbsPdf* constraintTerm = (RooAbsPdf*) constraints->find( constraintTermName.c_str() );
338  if( constraintTerm == NULL ) {
339  std::cout << "Error: Couldn't find constraint term: " << constraintTermName
340  << " for parameter: " << gamma_stat->GetName()
341  << std::endl;
342  throw std::runtime_error("Failed to find Gamma ConstraintTerm");
343  return -1;
344  }
345  */
346 
347  // Find the "data" of the poisson term
348  // This is the nominal value
349  bool FoundNomMean=false;
350  TIterator* iter_pois = constraintTerm->serverIterator(); //constraint_args
351  RooAbsArg* term_pois ;
352  while((term_pois=(RooAbsArg*)iter_pois->Next())) {
353  std::string serverName = term_pois->GetName();
354  //std::cout << "Checking Server: " << serverName << std::endl;
355  if( serverName.find("nom_")!=std::string::npos ) {
356  FoundNomMean = true;
357  pois_nom = (RooRealVar*) term_pois;
358  }
359  }
360  if( !FoundNomMean || !pois_nom ) {
361  std::cout << "Error: Did not find Nominal Pois Mean parameter in gamma constraint term PoissonMean: "
362  << constraintTerm->GetName() << std::endl;
363  throw std::runtime_error("Failed to find Nom Pois Mean");
364  }
365  else {
366  if(verbose) std::cout << "Found Poisson 'data' term: " << pois_nom->GetName() << std::endl;
367  }
368  delete iter_pois;
369 
370  // Taking the constraint term (a Poisson), find
371  // the "mean" which is the product: gamma*tau
372  // Then, from that mean, find tau
373  TIterator* iter_constr = constraintTerm->serverIterator(); //constraint_args
374  RooAbsArg* pois_mean_arg=NULL;
375  bool FoundPoissonMean = false;
376  while(( pois_mean_arg = (RooAbsArg*) iter_constr->Next() )) {
377  std::string serverName = pois_mean_arg->GetName();
378  if( pois_mean_arg->dependsOn( *gamma_stat ) ) {
379  FoundPoissonMean=true;
380  // pois_mean = (RooAbsReal*) pois_mean_arg;
381  break;
382  }
383  }
384  if( !FoundPoissonMean || !pois_mean_arg ) {
385  std::cout << "Error: Did not find PoissonMean parameter in gamma constraint term: "
386  << constraintTerm->GetName() << std::endl;
387  throw std::runtime_error("Failed to find PoissonMean");
388  return -1;
389  }
390  else {
391  if(verbose) std::cout << "Found Poisson 'mean' term: " << pois_mean_arg->GetName() << std::endl;
392  }
393  delete iter_constr;
394 
395 
396  TIterator* iter_product = pois_mean_arg->serverIterator(); //constraint_args
397  RooAbsArg* term_in_product ;
398  bool FoundTau=false;
399  while((term_in_product=(RooAbsArg*)iter_product->Next())) {
400  std::string serverName = term_in_product->GetName();
401  //std::cout << "Checking Server: " << serverName << std::endl;
402  if( serverName.find("_tau")!=std::string::npos ) {
403  FoundTau = true;
404  tau = (RooRealVar*) term_in_product;
405  }
406  }
407  if( !FoundTau || !tau ) {
408  std::cout << "Error: Did not find Tau parameter in gamma constraint term PoissonMean: "
409  << pois_mean_arg->GetName() << std::endl;
410  throw std::runtime_error("Failed to find Tau");
411  }
412  else {
413  if(verbose) std::cout << "Found Poisson 'tau' term: " << tau->GetName() << std::endl;
414  }
415  delete iter_product;
416 
417  return 0;
418 
419  }
420 
421 
422 
423 } // close RooStats namespace
424 } // close HistFactory namespace
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.
int getStatUncertaintyConstraintTerm(RooArgList *constraints, RooRealVar *gamma_stat, RooAbsReal *&pois_mean, RooRealVar *&tau)
const RooArgList & paramList() const
Definition: ParamHistFunc.h:39
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:744
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
virtual const RooArgSet * get() const
Definition: RooAbsData.h:77
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
bool getStatUncertaintyFromChannel(RooAbsPdf *channel, ParamHistFunc *&paramfunc, RooArgList *gammaList)
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4639
std::string channelNameFromPdf(RooAbsPdf *channelPdf)
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
#define NULL
Definition: RtypesCore.h:88
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:501
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:135
virtual Int_t numBins(const char *rangeName) const
Returm the number of fit bins ( = number of types )
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 tau
Definition: TRolke.cxx:630
TText * tt
Definition: textangle.C:16
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsArg.h:227
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state...
Definition: RooCatType.h:22
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
A doubly linked list.
Definition: TList.h:43
TIterator * typeIterator() const
Return iterator over all defined states.
virtual const Text_t * GetName() const
Returns name of object.
Definition: RooCatType.h:44
const RooAbsCategoryLValue & indexCat() const
RooAbsArg * first() const
const RooArgList & pdfList() const
Definition: RooProdPdf.h:69
RooArgSet * getComponents() const
Definition: RooAbsArg.cxx:688
bool verbose
void getDataValuesForObservables(std::map< std::string, std::vector< double > > &ChannelBinDataMap, RooAbsData *data, RooAbsPdf *simPdf)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
virtual const char * getLabel() const
Return label string of current state.
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
TIterator * serverIterator() const
Definition: RooAbsArg.h:112
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
The TH1 histogram class.
Definition: TH1.h:56
const char * GetName() const
Returns name of object.
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
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
virtual TList * split(const RooAbsCategory &splitCat, Bool_t createEmptyDataSets=kFALSE) const
Split dataset into subsets based on states of given splitCat in this dataset.
Bool_t contains(const RooAbsArg &var) const
virtual Int_t GetNbinsX() const
Definition: TH1.h:277
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
virtual void Print(Option_t *option="") const
Default print for collections, calls Print(option, 1).
RooAbsPdf * getSumPdfFromChannel(RooAbsPdf *channel)
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Int_t n
Definition: legend1.C:16
void FactorizeHistFactoryPdf(const RooArgSet &, RooAbsPdf &, RooArgList &, RooArgList &)