Logo ROOT  
Reference Guide
HistFactoryModelUtils.cxx
Go to the documentation of this file.
1/**
2 * \ingroup HistFactory
3 */
4
5// A set of utils for navegating HistFactory models
6#include <stdexcept>
7#include <typeinfo>
8
10#include "TIterator.h"
11#include "RooAbsArg.h"
12#include "RooAbsPdf.h"
13#include "RooArgSet.h"
14#include "RooArgList.h"
15#include "RooSimultaneous.h"
16#include "RooCategory.h"
17#include "RooRealVar.h"
18#include "RooProdPdf.h"
19#include "TH1.h"
20
23
24namespace RooStats{
25namespace HistFactory{
26
27
28 std::string channelNameFromPdf( RooAbsPdf* channelPdf ) {
29 std::string channelPdfName = channelPdf->GetName();
30 std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
31 return ChannelName;
32 }
33
35
36 bool verbose=false;
37
38 if(verbose) std::cout << "Getting the RooRealSumPdf for the channel: "
39 << sim_channel->GetName() << std::endl;
40
41 std::string channelPdfName = sim_channel->GetName();
42 std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
43
44 // Now, get the RooRealSumPdf
45 // ie the channel WITHOUT constraints
46 std::string realSumPdfName = ChannelName + "_model";
47
48 RooAbsPdf* sum_pdf = NULL;
49 TIterator* iter_sum_pdf = sim_channel->getComponents()->createIterator(); //serverIterator();
50 bool FoundSumPdf=false;
51 RooAbsArg* sum_pdf_arg=NULL;
52 while((sum_pdf_arg=(RooAbsArg*)iter_sum_pdf->Next())) {
53 std::string NodeClassName = sum_pdf_arg->ClassName();
54 if( NodeClassName == std::string("RooRealSumPdf") ) {
55 FoundSumPdf=true;
56 sum_pdf = (RooAbsPdf*) sum_pdf_arg;
57 break;
58 }
59 }
60 if( ! FoundSumPdf ) {
61 if(verbose) {
62 std::cout << "Failed to find RooRealSumPdf for channel: " << sim_channel->GetName() << std::endl;
63 sim_channel->getComponents()->Print("V");
64 }
65 sum_pdf=NULL;
66 //throw std::runtime_error("Failed to find RooRealSumPdf for channel");
67 }
68 else {
69 if(verbose) std::cout << "Found RooRealSumPdf: " << sum_pdf->GetName() << std::endl;
70 }
71 delete iter_sum_pdf;
72 iter_sum_pdf = NULL;
73
74 return sum_pdf;
75
76 }
77
78
79 void FactorizeHistFactoryPdf(const RooArgSet &observables, RooAbsPdf &pdf, RooArgList &obsTerms, RooArgList &constraints) {
80 // utility function to factorize constraint terms from a pdf
81 // (from G. Petrucciani)
82 const std::type_info & id = typeid(pdf);
83 if (id == typeid(RooProdPdf)) {
84 RooProdPdf *prod = dynamic_cast<RooProdPdf *>(&pdf);
85 RooArgList list(prod->pdfList());
86 for (int i = 0, n = list.getSize(); i < n; ++i) {
87 RooAbsPdf *pdfi = (RooAbsPdf *) list.at(i);
88 FactorizeHistFactoryPdf(observables, *pdfi, obsTerms, constraints);
89 }
90 } else if (id == typeid(RooSimultaneous) || id == typeid(HistFactorySimultaneous) ) { //|| id == typeid(RooSimultaneousOpt)) {
91 RooSimultaneous *sim = dynamic_cast<RooSimultaneous *>(&pdf);
93 for (int ic = 0, nc = cat->numBins((const char *)0); ic < nc; ++ic) {
94 cat->setBin(ic);
95 FactorizeHistFactoryPdf(observables, *sim->getPdf(cat->getCurrentLabel()), obsTerms, constraints);
96 }
97 delete cat;
98 } else if (pdf.dependsOn(observables)) {
99 if (!obsTerms.contains(pdf)) obsTerms.add(pdf);
100 } else {
101 if (!constraints.contains(pdf)) constraints.add(pdf);
102 }
103 }
104
105 /*
106 void getChannelsFromModel( RooAbsPdf* model, RooArgSet* channels, RooArgSet* channelsWithConstraints ) {
107
108 // Loop through the model
109 // Find all channels
110
111 std::string modelClassName = model->ClassName();
112
113 if( modelClassName == std::string("RooSimultaneous") || model->InheritsFrom("RooSimultaneous") ) {
114
115 TIterator* simServerItr = model->serverIterator();
116
117 // Loop through the child nodes of the sim pdf
118 // and find the channel nodes
119 RooAbsArg* sim_channel_arg = NULL;
120 while(( sim_channel = (RooAbsArg*) simServerItr->Next() )) {
121
122 RooAbsPdf* sim_channel = (RooAbsPdf*) sim_channel_arg;
123
124 // Ignore the Channel Cat
125 std::string channelPdfName = sim_channel->GetName();
126 std::string channelClassName = sim_channel->ClassName();
127 if( channelClassName == std::string("RooCategory") ) continue;
128
129 // If we got here, we found a channel.
130 // Format is model_<ChannelName>
131
132 std::string ChannelName = channelPdfName.substr(6, channelPdfName.size() );
133
134 // Now, get the RooRealSumPdf
135 RooAbsPdf* sum_pdf = getSumPdfFromChannel( sim_channel );
136
137
138 / *
139 // Now, get the RooRealSumPdf
140 // ie the channel WITHOUT constraints
141
142 std::string realSumPdfName = ChannelName + "_model";
143
144 RooAbsPdf* sum_pdf = NULL;
145 TIterator* iter_sum_pdf = sim_channel->getComponents()->createIterator(); //serverIterator();
146 bool FoundSumPdf=false;
147 RooAbsArg* sum_pdf_arg=NULL;
148 while((sum_pdf_arg=(RooAbsArg*)iter_sum_pdf->Next())) {
149
150 std::string NodeClassName = sum_pdf_arg->ClassName();
151 if( NodeClassName == std::string("RooRealSumPdf") ) {
152 FoundSumPdf=true;
153 sum_pdf = (RooAbsPdf*) sum_pdf_arg;
154 break;
155 }
156 }
157 if( ! FoundSumPdf ) {
158 std::cout << "Failed to find RooRealSumPdf for channel: " << sim_channel->GetName() << std::endl;
159 sim_channel->getComponents()->Print("V");
160 throw std::runtime_error("Failed to find RooRealSumPdf for channel");
161 }
162 delete iter_sum_pdf;
163 iter_sum_pdf = NULL;
164 * /
165
166 // Okay, now add to the arg sets
167 channels->add( *sum_pdf );
168 channelsWithConstraints->add( *sim_channel );
169
170 }
171
172 delete simServerItr;
173
174 }
175 else {
176 std::cout << "Model is not a RooSimultaneous or doesn't derive from one." << std::endl;
177 std::cout << "HistFactoryModelUtils isn't yet implemented for these pdf's" << std::endl;
178 }
179
180 }
181 */
182
183 bool getStatUncertaintyFromChannel( RooAbsPdf* channel, ParamHistFunc*& paramfunc, RooArgList* gammaList ) {
184
185 bool verbose=false;
186
187 // Find the servers of this channel
188 //TIterator* iter = channel->serverIterator();
189 TIterator* iter = channel->getComponents()->createIterator(); //serverIterator();
190 bool FoundParamHistFunc=false;
191 RooAbsArg* paramfunc_arg = NULL;
192 while(( paramfunc_arg = (RooAbsArg*) iter->Next() )) {
193 std::string NodeName = paramfunc_arg->GetName();
194 std::string NodeClassName = paramfunc_arg->ClassName();
195 if( NodeClassName != std::string("ParamHistFunc") ) continue;
196 if( NodeName.find("mc_stat_") != std::string::npos ) {
197 FoundParamHistFunc=true;
198 paramfunc = (ParamHistFunc*) paramfunc_arg;
199 break;
200 }
201 }
202 if( ! FoundParamHistFunc || !paramfunc ) {
203 if(verbose) std::cout << "Failed to find ParamHistFunc for channel: " << channel->GetName() << std::endl;
204 return false;
205 }
206
207 delete iter;
208 iter = NULL;
209
210 // Now, get the set of gamma's
211 gammaList = (RooArgList*) &( paramfunc->paramList());
212 if(verbose) gammaList->Print("V");
213
214 return true;
215
216 }
217
218
219 void getDataValuesForObservables( std::map< std::string, std::vector<double> >& ChannelBinDataMap,
220 RooAbsData* data, RooAbsPdf* pdf ) {
221
222 bool verbose=false;
223
224 //std::map< std::string, std::vector<int> ChannelBinDataMap;
225
226 RooSimultaneous* simPdf = (RooSimultaneous*) pdf;
227
228 // get category label
229 RooArgSet* allobs = (RooArgSet*) data->get();
230 TIterator* obsIter = allobs->createIterator();
231 RooCategory* cat = NULL;
232 RooAbsArg* temp = NULL;
233 while( (temp=(RooAbsArg*) obsIter->Next())) {
234 // use dynamic cast here instead
235 if( strcmp(temp->ClassName(),"RooCategory")==0){
236 cat = (RooCategory*) temp;
237 break;
238 }
239 }
240 if(verbose) {
241 if(!cat) std::cout <<"didn't find category"<< std::endl;
242 else std::cout <<"found category"<< std::endl;
243 }
244 delete obsIter;
245
246 // split dataset
247 TList* dataByCategory = data->split(*cat);
248 if(verbose) dataByCategory->Print();
249 // note :
250 // RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject("");
251
252 // loop over channels
253 RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
254 for (const auto& nameIdx : *channelCat) {
255
256 // Get pdf associated with state from simpdf
257 RooAbsPdf* pdftmp = simPdf->getPdf(nameIdx.first.c_str());
258
259 std::string ChannelName = pdftmp->GetName(); //tt->GetName();
260 if(verbose) std::cout << "Getting data for channel: " << ChannelName << std::endl;
261 ChannelBinDataMap[ ChannelName ] = std::vector<double>();
262
263 RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject(nameIdx.first.c_str());
264 if(verbose) dataForChan->Print();
265
266 // Generate observables defined by the pdf associated with this state
267 RooArgSet* obstmp = pdftmp->getObservables(*dataForChan->get()) ;
268 RooRealVar* obs = ((RooRealVar*)obstmp->first());
269 if(verbose) obs->Print();
270
271 //double expected = pdftmp->expectedEvents(*obstmp);
272
273 // set value to desired value (this is just an example)
274 // double obsVal = obs->getVal();
275 // set obs to desired value of observable
276 // obs->setVal( obsVal );
277 //double fracAtObsValue = pdftmp->getVal(*obstmp);
278
279 // get num events expected in bin for obsVal
280 // double nu = expected * fracAtObsValue;
281
282 // an easier way to get n
283 TH1* histForN = dataForChan->createHistogram("HhstForN",*obs);
284 for(int i=1; i<=histForN->GetNbinsX(); ++i){
285 double n = histForN->GetBinContent(i);
286 if(verbose) std::cout << "n" << i << " = " << n << std::endl;
287 ChannelBinDataMap[ ChannelName ].push_back( n );
288 }
289 delete histForN;
290
291 } // End Loop Over Categories
292
293 return;
294
295 }
296
297
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
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
Definition: ParamHistFunc.h:28
const RooArgList & paramList() const
Definition: ParamHistFunc.h:39
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:276
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:730
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:85
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:302
RooArgSet * getComponents() const
Definition: RooAbsArg.cxx:674
TIterator * serverIterator() const
Definition: RooAbsArg.h:141
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.
RooAbsArg * first() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
const char * GetName() const
Returns name of object.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
virtual TList * split(const RooAbsCategory &splitCat, Bool_t createEmptyDataSets=kFALSE) const
Split dataset into subsets based on states of given splitCat in this dataset.
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Definition: RooAbsData.cxx:628
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
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:74
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:23
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
const RooArgList & pdfList() const
Definition: RooProdPdf.h:67
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
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.
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
virtual void Print(Option_t *option="") const
Default print for collections, calls Print(option, 1).
The TH1 histogram class.
Definition: TH1.h:56
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
A doubly linked list.
Definition: TList.h:44
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:577
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
bool getStatUncertaintyFromChannel(RooAbsPdf *channel, ParamHistFunc *&paramfunc, RooArgList *gammaList)
void FactorizeHistFactoryPdf(const RooArgSet &, RooAbsPdf &, RooArgList &, RooArgList &)
void getDataValuesForObservables(std::map< std::string, std::vector< double > > &ChannelBinDataMap, RooAbsData *data, RooAbsPdf *simPdf)
RooAbsPdf * getSumPdfFromChannel(RooAbsPdf *channel)
int getStatUncertaintyConstraintTerm(RooArgList *constraints, RooRealVar *gamma_stat, RooAbsReal *&pois_mean, RooRealVar *&tau)
std::string channelNameFromPdf(RooAbsPdf *channelPdf)
const Int_t n
Definition: legend1.C:16
@ HistFactory
Definition: RooGlobalFunc.h:69
Namespace for the RooStats classes.
Definition: Asimov.h:19