Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 if (!cat) {
247 std::cerr <<"Category not found"<< std::endl;
248 return;
249 }
250
251 // split dataset
252 TList* dataByCategory = data->split(*cat);
253 if(verbose) dataByCategory->Print();
254 // note :
255 // RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject("");
256
257 // loop over channels
258 RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
259 for (const auto& nameIdx : *channelCat) {
260
261 // Get pdf associated with state from simpdf
262 RooAbsPdf* pdftmp = simPdf->getPdf(nameIdx.first.c_str());
263
264 std::string ChannelName = pdftmp->GetName(); //tt->GetName();
265 if(verbose) std::cout << "Getting data for channel: " << ChannelName << std::endl;
266 ChannelBinDataMap[ ChannelName ] = std::vector<double>();
267
268 RooAbsData* dataForChan = (RooAbsData*) dataByCategory->FindObject(nameIdx.first.c_str());
269 if(verbose) dataForChan->Print();
270
271 // Generate observables defined by the pdf associated with this state
272 RooArgSet* obstmp = pdftmp->getObservables(*dataForChan->get()) ;
273 RooRealVar* obs = ((RooRealVar*)obstmp->first());
274 if(verbose) obs->Print();
275
276 //double expected = pdftmp->expectedEvents(*obstmp);
277
278 // set value to desired value (this is just an example)
279 // double obsVal = obs->getVal();
280 // set obs to desired value of observable
281 // obs->setVal( obsVal );
282 //double fracAtObsValue = pdftmp->getVal(*obstmp);
283
284 // get num events expected in bin for obsVal
285 // double nu = expected * fracAtObsValue;
286
287 // an easier way to get n
288 TH1* histForN = dataForChan->createHistogram("HhstForN",*obs);
289 for(int i=1; i<=histForN->GetNbinsX(); ++i){
290 double n = histForN->GetBinContent(i);
291 if(verbose) std::cout << "n" << i << " = " << n << std::endl;
292 ChannelBinDataMap[ ChannelName ].push_back( n );
293 }
294 delete histForN;
295
296 } // End Loop Over Categories
297
298 dataByCategory->Delete();
299 delete dataByCategory;
300 }
301
302
304 RooAbsReal*& pois_nom, RooRealVar*& tau ) {
305 // Given a set of constraint terms,
306 // find the poisson constraint for the
307 // given gamma and return the mean
308 // as well as the 'tau' parameter
309
310 bool verbose=false;
311
312 // To get the constraint term, loop over all constraint terms
313 // and look for the gamma_stat name as well as '_constraint'
314 // std::string constraintTermName = std::string(gamma_stat->GetName()) + "_constraint";
315 TIterator* iter_list = constraints->createIterator();
316 RooAbsArg* term_constr=NULL;
317 bool FoundConstraintTerm=false;
318 RooAbsPdf* constraintTerm=NULL;
319 while((term_constr=(RooAbsArg*)iter_list->Next())) {
320 std::string TermName = term_constr->GetName();
321 // std::cout << "Checking if is a constraint term: " << TermName << std::endl;
322
323 //if( TermName.find(gamma_stat->GetName())!=string::npos ) {
324 if( term_constr->dependsOn( *gamma_stat) ) {
325 if( TermName.find("_constraint")!=std::string::npos ) {
326 FoundConstraintTerm=true;
327 constraintTerm = (RooAbsPdf*) term_constr;
328 break;
329 }
330 }
331 }
332 if( FoundConstraintTerm==false ) {
333 std::cout << "Error: Couldn't find constraint term for parameter: " << gamma_stat->GetName()
334 << " among constraints: " << constraints->GetName() << std::endl;
335 constraints->Print("V");
336 throw std::runtime_error("Failed to find Gamma ConstraintTerm");
337 return -1;
338 }
339 delete iter_list;
340
341 /*
342 RooAbsPdf* constraintTerm = (RooAbsPdf*) constraints->find( constraintTermName.c_str() );
343 if( constraintTerm == NULL ) {
344 std::cout << "Error: Couldn't find constraint term: " << constraintTermName
345 << " for parameter: " << gamma_stat->GetName()
346 << std::endl;
347 throw std::runtime_error("Failed to find Gamma ConstraintTerm");
348 return -1;
349 }
350 */
351
352 // Find the "data" of the poisson term
353 // This is the nominal value
354 bool FoundNomMean=false;
355 TIter iter_pois = constraintTerm->serverIterator(); //constraint_args
356 RooAbsArg* term_pois ;
357 while((term_pois=(RooAbsArg*)iter_pois.Next())) {
358 std::string serverName = term_pois->GetName();
359 //std::cout << "Checking Server: " << serverName << std::endl;
360 if( serverName.find("nom_")!=std::string::npos ) {
361 FoundNomMean = true;
362 pois_nom = (RooRealVar*) term_pois;
363 }
364 }
365 if( !FoundNomMean || !pois_nom ) {
366 std::cout << "Error: Did not find Nominal Pois Mean parameter in gamma constraint term PoissonMean: "
367 << constraintTerm->GetName() << std::endl;
368 throw std::runtime_error("Failed to find Nom Pois Mean");
369 }
370 else {
371 if(verbose) std::cout << "Found Poisson 'data' term: " << pois_nom->GetName() << std::endl;
372 }
373
374 // Taking the constraint term (a Poisson), find
375 // the "mean" which is the product: gamma*tau
376 // Then, from that mean, find tau
377 TIter iter_constr = constraintTerm->serverIterator(); //constraint_args
378 RooAbsArg* pois_mean_arg=NULL;
379 bool FoundPoissonMean = false;
380 while(( pois_mean_arg = (RooAbsArg*) iter_constr.Next() )) {
381 std::string serverName = pois_mean_arg->GetName();
382 if( pois_mean_arg->dependsOn( *gamma_stat ) ) {
383 FoundPoissonMean=true;
384 // pois_mean = (RooAbsReal*) pois_mean_arg;
385 break;
386 }
387 }
388 if( !FoundPoissonMean || !pois_mean_arg ) {
389 std::cout << "Error: Did not find PoissonMean parameter in gamma constraint term: "
390 << constraintTerm->GetName() << std::endl;
391 throw std::runtime_error("Failed to find PoissonMean");
392 return -1;
393 }
394 else {
395 if(verbose) std::cout << "Found Poisson 'mean' term: " << pois_mean_arg->GetName() << std::endl;
396 }
397
398 TIter iter_product = pois_mean_arg->serverIterator(); //constraint_args
399 RooAbsArg* term_in_product ;
400 bool FoundTau=false;
401 while((term_in_product=(RooAbsArg*)iter_product.Next())) {
402 std::string serverName = term_in_product->GetName();
403 //std::cout << "Checking Server: " << serverName << std::endl;
404 if( serverName.find("_tau")!=std::string::npos ) {
405 FoundTau = true;
406 tau = (RooRealVar*) term_in_product;
407 }
408 }
409 if( !FoundTau || !tau ) {
410 std::cout << "Error: Did not find Tau parameter in gamma constraint term PoissonMean: "
411 << pois_mean_arg->GetName() << std::endl;
412 throw std::runtime_error("Failed to find Tau");
413 }
414 else {
415 if(verbose) std::cout << "Found Poisson 'tau' term: " << tau->GetName() << std::endl;
416 }
417
418 return 0;
419
420 }
421
422
423
424} // close RooStats namespace
425} // close HistFactory namespace
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
const RooArgList & paramList() const
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:69
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition RooAbsArg.h:309
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:81
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:337
RooArgSet * getComponents() const
Create a RooArgSet with all components (branch nodes) of the expression tree headed by this object.
TIterator * serverIterator() const
Definition RooAbsArg.h:137
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=nullptr) 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:82
virtual const RooArgSet * get() const
Definition RooAbsData.h:128
virtual TList * split(const RooAbsCategory &splitCat, Bool_t createEmptyDataSets=kFALSE) const
Split dataset into subsets based on states of given splitCat in this dataset.
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition RooAbsData.h:250
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...
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:64
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:35
RooCategory is an object to represent discrete states.
Definition RooCategory.h:27
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
const RooArgList & pdfList() const
Definition RooProdPdf.h:69
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
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).
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Int_t GetNbinsX() const
Definition TH1.h:296
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4994
TObject * Next()
Iterator abstract base class.
Definition TIterator.h:30
virtual TObject * Next()=0
A doubly linked list.
Definition TList.h:38
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition TList.cxx:578
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:470
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:200
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
Namespace for the RooStats classes.
Definition Asimov.h:19