Logo ROOT  
Reference Guide
RooAddModel.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17//////////////////////////////////////////////////////////////////////////////
18///
19/// RooAddModel is an efficient implementation of a sum of PDFs of the form
20/// \f[
21/// c_1*\mathrm{PDF}_1 + c_2*\mathrm{PDF}_2 + ... c_n*\mathrm{PDF}_n
22/// \]f
23/// or
24/// \f[
25/// c_1*\mathrm{PDF}_1 + c_2*\mathrm{PDF}_2 + ... (1-\sum(c_1, \ldots, c_{n-1}))*\mathrm{PDF}_n
26/// \f]
27/// The first form is for extended likelihood fits, where the
28/// expected number of events is \f$ \sum_i c_i \f$. The coefficients \f$ c_i \f$
29/// can either be explicitly provided, or, if all components support
30/// extended likelihood fits, they can be calculated the contribution
31/// of each PDF to the total number of expected events.
32///
33/// In the second form, the sum of the coefficients is enforced to be one,
34/// and the coefficient of the last PDF is calculated from that condition.
35///
36/// RooAddPdf relies on each component PDF to be normalized and will perform
37/// no normalization other than calculating the proper last coefficient \f$ c_n \f$, if requested.
38/// An (enforced) condition for this assumption is that each \f$ \mathrm{PDF}_i \f$ is independent
39/// of each coefficient i.
40///
41///
42
43#include "RooAddModel.h"
44
45#include "RooFit.h"
46#include "RooMsgService.h"
47#include "RooDataSet.h"
48#include "RooRealProxy.h"
49#include "RooPlot.h"
50#include "RooRealVar.h"
51#include "RooAddGenContext.h"
52#include "RooRealConstant.h"
53#include "RooNameReg.h"
54#include "RooRealIntegral.h"
55
56using namespace std;
57
59;
60
61
62////////////////////////////////////////////////////////////////////////////////
63
65 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
66 _refCoefRangeName(0),
67 _projectCoefs(false),
68 _codeReg(10),
69 _snormList(0),
70 _haveLastCoef(false),
71 _allExtendable(false)
72{
73 _coefCache = new Double_t[10] ;
75}
76
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// Generic constructor from list of PDFs and list of coefficients.
81/// Each pdf list element (i) is paired with coefficient list element (i).
82/// The number of coefficients must be either equal to the number of PDFs,
83/// in which case extended MLL fitting is enabled, or be one less.
84///
85/// All PDFs must inherit from RooAbsPdf. All coefficients must inherit from RooAbsReal
86
87RooAddModel::RooAddModel(const char *name, const char *title, const RooArgList& inPdfList, const RooArgList& inCoefList, Bool_t ownPdfList) :
88 RooResolutionModel(name,title,((RooResolutionModel*)inPdfList.at(0))->convVar()),
89 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
90 _refCoefRangeName(0),
91 _projectCoefs(kFALSE),
92 _projCacheMgr(this,10),
93 _intCacheMgr(this,10),
94 _codeReg(10),
95 _pdfList("!pdfs","List of PDFs",this),
96 _coefList("!coefficients","List of coefficients",this),
97 _haveLastCoef(kFALSE),
98 _allExtendable(kFALSE)
99{
100 if (inPdfList.getSize()>inCoefList.getSize()+1) {
101 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName()
102 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
103 assert(0) ;
104 }
105
106 // Constructor with N PDFs and N or N-1 coefs
107 TIterator* pdfIter = inPdfList.createIterator() ;
108 TIterator* coefIter = inCoefList.createIterator() ;
109 RooAbsPdf* pdf ;
110 RooAbsReal* coef ;
111
112 while((coef = (RooAbsPdf*)coefIter->Next())) {
113 pdf = (RooAbsPdf*) pdfIter->Next() ;
114 if (!pdf) {
115 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName()
116 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
117 assert(0) ;
118 }
119 if (!dynamic_cast<RooAbsReal*>(coef)) {
120 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal, ignored" << endl ;
121 continue ;
122 }
123 if (!dynamic_cast<RooAbsReal*>(pdf)) {
124 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") pdf " << pdf->GetName() << " is not of type RooAbsPdf, ignored" << endl ;
125 continue ;
126 }
127 _pdfList.add(*pdf) ;
128 _coefList.add(*coef) ;
129 }
130
131 pdf = (RooAbsPdf*) pdfIter->Next() ;
132 if (pdf) {
133 if (!dynamic_cast<RooAbsReal*>(pdf)) {
134 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") last pdf " << coef->GetName() << " is not of type RooAbsPdf, fatal error" << endl ;
135 assert(0) ;
136 }
137 _pdfList.add(*pdf) ;
138 } else {
140 }
141
142 delete pdfIter ;
143 delete coefIter ;
144
147
148 if (ownPdfList) {
150 }
151
152}
153
154
155
156////////////////////////////////////////////////////////////////////////////////
157/// Copy constructor
158
159RooAddModel::RooAddModel(const RooAddModel& other, const char* name) :
161 _refCoefNorm("!refCoefNorm",this,other._refCoefNorm),
162 _refCoefRangeName((TNamed*)other._refCoefRangeName),
163 _projectCoefs(other._projectCoefs),
164 _projCacheMgr(other._projCacheMgr,this),
165 _intCacheMgr(other._intCacheMgr,this),
166 _codeReg(other._codeReg),
167 _pdfList("!pdfs",this,other._pdfList),
168 _coefList("!coefficients",this,other._coefList),
169 _haveLastCoef(other._haveLastCoef),
170 _allExtendable(other._allExtendable)
171{
174}
175
176
177
178////////////////////////////////////////////////////////////////////////////////
179/// Destructor
180
182{
183 if (_coefCache) delete[] _coefCache ;
184}
185
186
187
188////////////////////////////////////////////////////////////////////////////////
189/// By default the interpretation of the fraction coefficients is
190/// performed in the contextual choice of observables. This makes the
191/// shape of the p.d.f explicitly dependent on the choice of
192/// observables. This method instructs RooAddPdf to freeze the
193/// interpretation of the coefficients to be done in the given set of
194/// observables. If frozen, fractions are automatically transformed
195/// from the reference normalization set to the contextual normalization
196/// set by ratios of integrals
197
199{
200 if (refCoefNorm.getSize()==0) {
202 return ;
203 }
205
207 _refCoefNorm.add(refCoefNorm) ;
208
210}
211
212
213
214////////////////////////////////////////////////////////////////////////////////
215/// By default the interpretation of the fraction coefficients is
216/// performed in the default range. This make the shape of a RooAddPdf
217/// explicitly dependent on the range of the observables. To allow
218/// a range independent definition of the fraction this function
219/// instructs RooAddPdf to freeze its interpretation in the given
220/// named range. If the current normalization range is different
221/// from the reference range, the appropriate fraction coefficients
222/// are automically calculation from the reference fractions using
223/// ratios if integrals
224
225void RooAddModel::fixCoefRange(const char* rangeName)
226{
229}
230
231
232
233////////////////////////////////////////////////////////////////////////////////
234/// Instantiate a clone of this resolution model representing a convolution with given
235/// basis function. The owners object name is incorporated in the clones name
236/// to avoid multiple convolution objects with the same name in complex PDF structures.
237///
238/// RooAddModel will clone all the component models to create a composite convolution object
239
241{
242 // Check that primary variable of basis functions is our convolution variable
243 if (inBasis->getParameter(0) != x.absArg()) {
244 coutE(InputArguments) << "RooAddModel::convolution(" << GetName()
245 << ") convolution parameter of basis function and PDF don't match" << endl ;
246 ccoutE(InputArguments) << "basis->findServer(0) = " << inBasis->findServer(0) << " " << inBasis->findServer(0)->GetName() << endl ;
247 ccoutE(InputArguments) << "x.absArg() = " << x.absArg() << " " << x.absArg()->GetName() << endl ;
248 inBasis->Print("v") ;
249 return 0 ;
250 }
251
252 TString newName(GetName()) ;
253 newName.Append("_conv_") ;
254 newName.Append(inBasis->GetName()) ;
255 newName.Append("_[") ;
256 newName.Append(owner->GetName()) ;
257 newName.Append("]") ;
258
259 TString newTitle(GetTitle()) ;
260 newTitle.Append(" convoluted with basis function ") ;
261 newTitle.Append(inBasis->GetName()) ;
262
263 RooArgList modelList ;
264 for (auto obj : _pdfList) {
265 auto model = static_cast<RooResolutionModel*>(obj);
266 // Create component convolution
267 RooResolutionModel* conv = model->convolution(inBasis,owner) ;
268 modelList.add(*conv) ;
269 }
270
271 RooArgList theCoefList ;
272 for (auto coef : _coefList) {
273 theCoefList.add(*coef) ;
274 }
275
276 RooAddModel* convSum = new RooAddModel(newName,newTitle,modelList,theCoefList,kTRUE) ;
277 for (std::set<std::string>::const_iterator attrIt = _boolAttrib.begin();
278 attrIt != _boolAttrib.end(); ++attrIt) {
279 convSum->setAttribute((*attrIt).c_str()) ;
280 }
281 for (std::map<std::string,std::string>::const_iterator attrIt = _stringAttrib.begin();
282 attrIt != _stringAttrib.end(); ++attrIt) {
283 convSum->setStringAttribute((attrIt->first).c_str(), (attrIt->second).c_str()) ;
284 }
285 convSum->changeBasis(inBasis) ;
286 return convSum ;
287}
288
289
290
291////////////////////////////////////////////////////////////////////////////////
292/// Return code for basis function representing by 'name' string.
293/// The basis code of the first component model will be returned,
294/// if the basis is supported by all components. Otherwise 0
295/// is returned
296
298{
299 Bool_t first(kTRUE), code(0) ;
300 for (auto obj : _pdfList) {
301 auto model = static_cast<RooResolutionModel*>(obj);
302 Int_t subCode = model->basisCode(name) ;
303 if (first) {
304 code = subCode ;
305 first = kFALSE ;
306 } else if (subCode==0) {
307 code = 0 ;
308 }
309 }
310
311 return code ;
312}
313
314
315
316////////////////////////////////////////////////////////////////////////////////
317/// Retrieve cache element with for calculation of p.d.f value with normalization set nset and integrated over iset
318/// in range 'rangeName'. If cache element does not exist, create and fill it on the fly. The cache contains
319/// suplemental normalization terms (in case not all added p.d.f.s have the same observables), projection
320/// integrals to calculated transformed fraction coefficients when a frozen reference frame is provided
321/// and projection integrals for similar transformations when a frozen reference range is provided.
322
323RooAddModel::CacheElem* RooAddModel::getProjCache(const RooArgSet* nset, const RooArgSet* iset, const char* rangeName) const
324{
325 // Check if cache already exists
326 CacheElem* cache = (CacheElem*) _projCacheMgr.getObj(nset,iset,0,RooNameReg::ptr(rangeName)) ;
327 if (cache) {
328 return cache ;
329 }
330
331 //Create new cache
332 cache = new CacheElem ;
333
334 // *** PART 1 : Create supplemental normalization list ***
335
336 // Retrieve the combined set of dependents of this PDF ;
337 RooArgSet *fullDepList = getObservables(nset) ;
338 if (iset) {
339 fullDepList->remove(*iset,kTRUE,kTRUE) ;
340 }
341
342 // Fill with dummy unit RRVs for now
343 for (unsigned int i = 0; i < _pdfList.size(); ++i) {
344 auto pdf = static_cast<RooAbsPdf*>(&_pdfList[i]);
345 auto coef = i < _coefList.size() ? &_coefList[i] : nullptr;
346
347 // Start with full list of dependents
348 RooArgSet supNSet(*fullDepList) ;
349
350 // Remove PDF dependents
351 RooArgSet* pdfDeps = pdf->getObservables(nset) ;
352 if (pdfDeps) {
353 supNSet.remove(*pdfDeps,kTRUE,kTRUE) ;
354 delete pdfDeps ;
355 }
356
357 // Remove coef dependents
358 if (coef) {
359 RooArgSet* coefDeps = coef->getObservables(nset);
360 supNSet.remove(*coefDeps,kTRUE,kTRUE) ;
361 delete coefDeps ;
362 }
363
364 RooAbsReal* snorm ;
365 TString name(GetName()) ;
366 name.Append("_") ;
367 name.Append(pdf->GetName()) ;
368 name.Append("_SupNorm") ;
369 if (supNSet.getSize()>0) {
370 snorm = new RooRealIntegral(name,"Supplemental normalization integral",RooRealConstant::value(1.0),supNSet) ;
371 } else {
372 snorm = new RooRealVar(name,"Unit Supplemental normalization integral",1.0) ;
373 }
374 cache->_suppNormList.addOwned(*snorm) ;
375 }
376
377 delete fullDepList ;
378
379 if (_verboseEval>1) {
380 cxcoutD(Caching) << "RooAddModel::syncSuppNormList(" << GetName() << ") synching supplemental normalization list for norm" << (nset?*nset:RooArgSet()) << endl ;
381 if (dologD(Caching)) {
382 cache->_suppNormList.Print("v") ;
383 }
384 }
385
386
387 // *** PART 2 : Create projection coefficients ***
388
389 // If no projections required stop here
390 if (!_projectCoefs || _basis!=0 ) {
391 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
392 return cache ;
393 }
394
395
396 // Reduce iset/nset to actual dependents of this PDF
397 RooArgSet* nset2 = nset ? getObservables(nset) : new RooArgSet() ;
398
399 // Check if requested transformation is not identity
400 if (!nset2->equals(_refCoefNorm) || _refCoefRangeName !=0 || rangeName !=0 ) {
401
402 coutI(Caching) << "RooAddModel::syncCoefProjList(" << GetName() << ") creating coefficient projection integrals" << endl ;
403 ccoutI(Caching) << " from current normalization: " ; nset2->Print("1") ;
404 ccoutI(Caching) << " with current range: " << (rangeName?rangeName:"<none>") << endl ;
405 ccoutI(Caching) << " to reference normalization: " ; _refCoefNorm.Print("1") ;
406 ccoutI(Caching) << " with reference range: " << (_refCoefRangeName?RooNameReg::str(_refCoefRangeName):"<none>") << endl ;
407
408 // Recalculate projection integrals of PDFs
409 for (const auto obj : _pdfList) {
410 const auto thePdf = static_cast<RooAbsPdf*>(obj);
411
412 // Calculate projection integral
413 RooAbsReal* pdfProj ;
414 if (!nset2->equals(_refCoefNorm)) {
415 pdfProj = thePdf->createIntegral(*nset2,_refCoefNorm) ;
416 pdfProj->setOperMode(operMode()) ;
417 } else {
418 TString name(GetName()) ;
419 name.Append("_") ;
420 name.Append(thePdf->GetName()) ;
421 name.Append("_ProjectNorm") ;
422 pdfProj = new RooRealVar(name,"Unit Projection normalization integral",1.0) ;
423 }
424
425 cache->_projList.addOwned(*pdfProj) ;
426
427 // Calculation optional supplemental normalization term
428 RooArgSet supNormSet(_refCoefNorm) ;
429 RooArgSet* deps = thePdf->getParameters(RooArgSet()) ;
430 supNormSet.remove(*deps,kTRUE,kTRUE) ;
431 delete deps ;
432
433 RooAbsReal* snorm ;
434 TString name(GetName()) ;
435 name.Append("_") ;
436 name.Append(thePdf->GetName()) ;
437 name.Append("_ProjSupNorm") ;
438 if (supNormSet.getSize()>0) {
439 snorm = new RooRealIntegral(name,"Projection Supplemental normalization integral",
440 RooRealConstant::value(1.0),supNormSet) ;
441 } else {
442 snorm = new RooRealVar(name,"Unit Projection Supplemental normalization integral",1.0) ;
443 }
444 cache->_suppProjList.addOwned(*snorm) ;
445
446 // Calculate reference range adjusted projection integral
447 RooAbsReal* rangeProj1 ;
450 rangeProj1->setOperMode(operMode()) ;
451 } else {
452 TString theName(GetName()) ;
453 theName.Append("_") ;
454 theName.Append(thePdf->GetName()) ;
455 theName.Append("_RangeNorm1") ;
456 rangeProj1 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
457 }
458 cache->_refRangeProjList.addOwned(*rangeProj1) ;
459
460
461 // Calculate range adjusted projection integral
462 RooAbsReal* rangeProj2 ;
463 if (rangeName && _refCoefNorm.getSize()>0) {
464 rangeProj2 = thePdf->createIntegral(_refCoefNorm,_refCoefNorm,rangeName) ;
465 rangeProj2->setOperMode(operMode()) ;
466 } else {
467 TString theName(GetName()) ;
468 theName.Append("_") ;
469 theName.Append(thePdf->GetName()) ;
470 theName.Append("_RangeNorm2") ;
471 rangeProj2 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
472 }
473 cache->_rangeProjList.addOwned(*rangeProj2) ;
474
475 }
476
477 }
478
479 delete nset2 ;
480
481 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
482
483 return cache ;
484}
485
486
487
488////////////////////////////////////////////////////////////////////////////////
489/// Update the coefficient values in the given cache element: calculate new remainder
490/// fraction, normalize fractions obtained from extended ML terms to unity and
491/// multiply these the various range and dimensional corrections needed in the
492/// current use context
493
495{
496 // cxcoutD(ChangeTracking) << "RooAddModel::updateCoefficients(" << GetName() << ") update coefficients" << endl ;
497
498 Int_t i ;
499
500 // Straight coefficients
501 if (_allExtendable) {
502
503 // coef[i] = expectedEvents[i] / SUM(expectedEvents)
504 Double_t coefSum(0) ;
505 for (i=0 ; i<_pdfList.getSize() ; i++) {
507 coefSum += _coefCache[i] ;
508 }
509 if (coefSum==0.) {
510 coutW(Eval) << "RooAddModel::updateCoefCache(" << GetName() << ") WARNING: total number of expected events is 0" << endl ;
511 } else {
512 for (i=0 ; i<_pdfList.getSize() ; i++) {
513 _coefCache[i] /= coefSum ;
514 }
515 }
516
517 } else {
518 if (_haveLastCoef) {
519
520 // coef[i] = coef[i] / SUM(coef)
521 Double_t coefSum(0) ;
522 for (i=0 ; i<_coefList.getSize() ; i++) {
523 _coefCache[i] = ((RooAbsPdf*)_coefList.at(i))->getVal(nset) ;
524 coefSum += _coefCache[i] ;
525 }
526 for (i=0 ; i<_coefList.getSize() ; i++) {
527 _coefCache[i] /= coefSum ;
528 }
529 } else {
530
531 // coef[i] = coef[i] ; coef[n] = 1-SUM(coef[0...n-1])
532 Double_t lastCoef(1) ;
533 for (i=0 ; i<_coefList.getSize() ; i++) {
534 _coefCache[i] = ((RooAbsPdf*)_coefList.at(i))->getVal(nset) ;
535 cxcoutD(Caching) << "SYNC: orig coef[" << i << "] = " << _coefCache[i] << endl ;
536 lastCoef -= _coefCache[i] ;
537 }
538 _coefCache[_coefList.getSize()] = lastCoef ;
539 cxcoutD(Caching) << "SYNC: orig coef[" << _coefList.getSize() << "] = " << _coefCache[_coefList.getSize()] << endl ;
540
541
542 // Warn about coefficient degeneration
543 if ((lastCoef<-1e-05 || (lastCoef-1)>1e-5) && _coefErrCount-->0) {
544 coutW(Eval) << "RooAddModel::updateCoefCache(" << GetName()
545 << " WARNING: sum of PDF coefficients not in range [0-1], value="
546 << 1-lastCoef << endl ;
547 if (_coefErrCount==0) {
548 coutW(Eval) << " (no more will be printed)" << endl ;
549 }
550 }
551 }
552 }
553
554
555
556 // Stop here if not projection is required or needed
557 if ((!_projectCoefs) || cache._projList.getSize()==0) {
558 // cout << "SYNC no projection required rangeName = " << (rangeName?rangeName:"<none>") << endl ;
559 return ;
560 }
561
562 // Adjust coefficients for given projection
563 Double_t coefSum(0) ;
564 for (i=0 ; i<_pdfList.getSize() ; i++) {
565 GlobalSelectComponentRAII compRAII(true);
566
567 RooAbsReal* pp = ((RooAbsReal*)cache._projList.at(i)) ;
568 RooAbsReal* sn = ((RooAbsReal*)cache._suppProjList.at(i)) ;
569 RooAbsReal* r1 = ((RooAbsReal*)cache._refRangeProjList.at(i)) ;
570 RooAbsReal* r2 = ((RooAbsReal*)cache._rangeProjList.at(i)) ;
571
572 if (dologD(Eval)) {
573 cxcoutD(Eval) << "pp = " << pp->GetName() << endl
574 << "sn = " << sn->GetName() << endl
575 << "r1 = " << r1->GetName() << endl
576 << "r2 = " << r2->GetName() << endl ;
579 }
580
581 Double_t proj = pp->getVal()/sn->getVal()*(r2->getVal()/r1->getVal()) ;
582
583 _coefCache[i] *= proj ;
584 coefSum += _coefCache[i] ;
585 }
586 for (i=0 ; i<_pdfList.getSize() ; i++) {
587 _coefCache[i] /= coefSum ;
588// cout << "POST-SYNC coef[" << i << "] = " << _coefCache[i] << endl ;
589 }
590
591
592
593}
594
595
596
597////////////////////////////////////////////////////////////////////////////////
598/// Calculate the current value
599
601{
602 const RooArgSet* nset = _normSet ;
603 CacheElem* cache = getProjCache(nset) ;
604
605 updateCoefficients(*cache,nset) ;
606
607
608 // Do running sum of coef/pdf pairs, calculate lastCoef.
609 Double_t snormVal ;
610 Double_t value(0) ;
611 Int_t i(0) ;
612 for (auto obj : _pdfList) {
613 auto pdf = static_cast<RooAbsPdf*>(obj);
614
615 if (_coefCache[i]!=0.) {
616 snormVal = nset ? ((RooAbsReal*)cache->_suppNormList.at(i))->getVal() : 1.0 ;
617 Double_t pdfVal = pdf->getVal(nset) ;
618 // Double_t pdfNorm = pdf->getNorm(nset) ;
619 if (pdf->isSelectedComp()) {
620 value += pdfVal*_coefCache[i]/snormVal ;
621 cxcoutD(Eval) << "RooAddModel::evaluate(" << GetName() << ") value += ["
622 << pdf->GetName() << "] " << pdfVal << " * " << _coefCache[i] << " / " << snormVal << endl ;
623 }
624 }
625 i++ ;
626 }
627
628 return value ;
629}
630
631
632
633////////////////////////////////////////////////////////////////////////////////
634/// Reset error counter to given value, limiting the number
635/// of future error messages for this pdf to 'resetValue'
636
638{
640 _coefErrCount = resetValue ;
641}
642
643
644
645////////////////////////////////////////////////////////////////////////////////
646/// Check if PDF is valid for given normalization set.
647/// Coeffient and PDF must be non-overlapping, but pdf-coefficient
648/// pairs may overlap each other
649
651{
652 Bool_t ret(kFALSE) ;
653
654 for (unsigned int i = 0; i < _coefList.size(); ++i) {
655 auto pdf = &_pdfList[i];
656 auto coef = &_coefList[i];
657
658 if (pdf->observableOverlaps(nset,*coef)) {
659 coutE(InputArguments) << "RooAddModel::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName()
660 << " and PDF " << pdf->GetName() << " have one or more dependents in common" << endl ;
661 ret = kTRUE ;
662 }
663 }
664
665 return ret ;
666}
667
668
669
670////////////////////////////////////////////////////////////////////////////////
671
673 const RooArgSet* normSet, const char* rangeName) const
674{
675 if (_forceNumInt) return 0 ;
676
677 // Declare that we can analytically integrate all requested observables
678 analVars.add(allVars) ;
679
680 // Retrieve (or create) the required component integral list
681 Int_t code ;
682 RooArgList *cilist ;
683 getCompIntList(normSet,&allVars,cilist,code,rangeName) ;
684
685 return code+1 ;
686
687}
688
689
690
691////////////////////////////////////////////////////////////////////////////////
692/// Check if this configuration was created before
693
694void RooAddModel::getCompIntList(const RooArgSet* nset, const RooArgSet* iset, pRooArgList& compIntList, Int_t& code, const char* isetRangeName) const
695{
696 Int_t sterileIdx(-1) ;
697
698 IntCacheElem* cache = (IntCacheElem*) _intCacheMgr.getObj(nset,iset,&sterileIdx,RooNameReg::ptr(isetRangeName)) ;
699 if (cache) {
700 code = _intCacheMgr.lastIndex() ;
701 compIntList = &cache->_intList ;
702
703 return ;
704 }
705
706 // Create containers for partial integral components to be generated
707 cache = new IntCacheElem ;
708
709 // Fill Cache
710 for (auto obj : _pdfList) {
711 auto model = static_cast<RooResolutionModel*>(obj);
712
713 RooAbsReal* intPdf = model->createIntegral(*iset,nset,0,isetRangeName) ;
714 cache->_intList.addOwned(*intPdf) ;
715 }
716
717 // Store the partial integral list and return the assigned code ;
718 code = _intCacheMgr.setObj(nset,iset,(RooAbsCacheElement*)cache,RooNameReg::ptr(isetRangeName)) ;
719
720 // Fill references to be returned
721 compIntList = &cache->_intList ;
722}
723
724
725
726////////////////////////////////////////////////////////////////////////////////
727/// Return analytical integral defined by given scenario code
728
729Double_t RooAddModel::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
730{
731 // No integration scenario
732 if (code==0) {
733 return getVal(normSet) ;
734 }
735
736 // Partial integration scenarios
738
739 RooArgList* compIntList ;
740
741 // If cache has been sterilized, revive this slot
742 if (cache==0) {
744 RooArgSet* nset = _intCacheMgr.nameSet1ByIndex(code-1)->select(*vars) ;
745 RooArgSet* iset = _intCacheMgr.nameSet2ByIndex(code-1)->select(*vars) ;
746
747 Int_t code2(-1) ;
748 getCompIntList(nset,iset,compIntList,code2,rangeName) ;
749
750 delete vars ;
751 delete nset ;
752 delete iset ;
753 } else {
754
755 compIntList = &cache->_intList ;
756
757 }
758
759 // Calculate the current value
760 const RooArgSet* nset = _normSet ;
761 CacheElem* pcache = getProjCache(nset) ;
762
763 updateCoefficients(*pcache,nset) ;
764
765 // Do running sum of coef/pdf pairs, calculate lastCoef.
766 Double_t snormVal ;
767 Double_t value(0) ;
768 Int_t i(0) ;
769 for (const auto obj : *compIntList) {
770 auto pdfInt = static_cast<const RooAbsReal*>(obj);
771 if (_coefCache[i]!=0.) {
772 snormVal = nset ? ((RooAbsReal*)pcache->_suppNormList.at(i))->getVal() : 1.0 ;
773 Double_t intVal = pdfInt->getVal(nset) ;
774 value += intVal*_coefCache[i]/snormVal ;
775 cxcoutD(Eval) << "RooAddModel::evaluate(" << GetName() << ") value += ["
776 << pdfInt->GetName() << "] " << intVal << " * " << _coefCache[i] << " / " << snormVal << endl ;
777 }
778 i++ ;
779 }
780
781 return value ;
782
783}
784
785
786
787////////////////////////////////////////////////////////////////////////////////
788/// Return the number of expected events, which is either the sum of all coefficients
789/// or the sum of the components extended terms
790
792{
793 Double_t expectedTotal(0.0);
794
795 if (_allExtendable) {
796
797 // Sum of the extended terms
798 for (auto obj : _pdfList) {
799 auto pdf = static_cast<RooAbsPdf*>(obj);
800 expectedTotal += pdf->expectedEvents(nset) ;
801 }
802
803 } else {
804
805 // Sum the coefficients
806 for (const auto obj : _coefList) {
807 auto coef = static_cast<RooAbsReal*>(obj);
808 expectedTotal += coef->getVal() ;
809 }
810 }
811
812 return expectedTotal;
813}
814
815
816
817////////////////////////////////////////////////////////////////////////////////
818/// Interface function used by test statistics to freeze choice of observables
819/// for interpretation of fraction coefficients
820
822{
823 if (!force && _refCoefNorm.getSize()!=0) {
824 return ;
825 }
826
827 if (!depSet) {
829 return ;
830 }
831
832 RooArgSet* myDepSet = getObservables(depSet) ;
833 fixCoefNormalization(*myDepSet) ;
834 delete myDepSet ;
835}
836
837
838
839////////////////////////////////////////////////////////////////////////////////
840/// Interface function used by test statistics to freeze choice of range
841/// for interpretation of fraction coefficients
842
843void RooAddModel::selectNormalizationRange(const char* rangeName, Bool_t force)
844{
845 if (!force && _refCoefRangeName) {
846 return ;
847 }
848
849 fixCoefRange(rangeName) ;
850}
851
852
853
854////////////////////////////////////////////////////////////////////////////////
855/// Return specialized context to efficiently generate toy events from RooAddPdfs
856
858 const RooArgSet* auxProto, Bool_t verbose) const
859{
860 return new RooAddGenContext(*this,vars,prototype,auxProto,verbose) ;
861}
862
863
864
865////////////////////////////////////////////////////////////////////////////////
866/// Direct generation is safe if all components say so
867
869{
870 for (auto obj : _pdfList) {
871 auto pdf = static_cast<RooAbsPdf*>(obj);
872
873 if (!pdf->isDirectGenSafe(arg)) {
874 return kFALSE ;
875 }
876 }
877 return kTRUE ;
878}
879
880
881
882////////////////////////////////////////////////////////////////////////////////
883/// Return pseud-code that indicates if all components can do internal generation (1) or not (0)
884
885Int_t RooAddModel::getGenerator(const RooArgSet& directVars, RooArgSet &/*generateVars*/, Bool_t /*staticInitOK*/) const
886{
887 for (auto obj : _pdfList) {
888 auto pdf = static_cast<RooAbsPdf*>(obj);
889
890 RooArgSet tmp ;
891 if (pdf->getGenerator(directVars,tmp)==0) {
892 return 0 ;
893 }
894 }
895 return 1 ;
896}
897
898
899
900
901////////////////////////////////////////////////////////////////////////////////
902/// This function should never be called as RooAddModel implements a custom generator context
903
905{
906 assert(0) ;
907}
908
909
910
911
912////////////////////////////////////////////////////////////////////////////////
913/// List all RooAbsArg derived contents in this cache element
914
916{
917 RooArgList allNodes;
918 allNodes.add(_projList) ;
919 allNodes.add(_suppProjList) ;
920 allNodes.add(_refRangeProjList) ;
921 allNodes.add(_rangeProjList) ;
922
923 return allNodes ;
924}
925
926
927
928////////////////////////////////////////////////////////////////////////////////
929/// List all RooAbsArg derived contents in this cache element
930
932{
933 RooArgList allNodes(_intList) ;
934 return allNodes ;
935}
936
937
938////////////////////////////////////////////////////////////////////////////////
939/// Customized printing of arguments of a RooAddModel to more intuitively reflect the contents of the
940/// product operator construction
941
942void RooAddModel::printMetaArgs(ostream& os) const
943{
945
946 os << "(" ;
947 for (unsigned int i=0; i < _coefList.size(); ++i) {
948 auto coef = &_coefList[i];
949 auto pdf = &_pdfList[i];
950 if (!first) {
951 os << " + " ;
952 } else {
953 first = kFALSE ;
954 }
955 os << coef->GetName() << " * " << pdf->GetName() ;
956 }
957 if (_pdfList.size() > _coefList.size()) {
958 os << " + [%] * " << _pdfList[_pdfList.size()-1].GetName() ;
959 }
960 os << ") " ;
961}
962
#define e(i)
Definition: RSha256.hxx:103
#define coutI(a)
Definition: RooMsgService.h:30
#define ccoutE(a)
Definition: RooMsgService.h:41
#define cxcoutD(a)
Definition: RooMsgService.h:81
#define coutW(a)
Definition: RooMsgService.h:32
#define dologD(a)
Definition: RooMsgService.h:65
#define coutE(a)
Definition: RooMsgService.h:33
#define ccoutI(a)
Definition: RooMsgService.h:38
#define ccoutD(a)
Definition: RooMsgService.h:37
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
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
void printCompactTree(const char *indent="", const char *fileName=0, const char *namePat=0, RooAbsArg *client=0)
Print tree structure of expression tree on stdout, or to file if filename is specified.
Definition: RooAbsArg.cxx:1745
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
Definition: RooAbsArg.cxx:289
friend class RooArgSet
Definition: RooAbsArg.h:572
std::set< std::string > _boolAttrib
Definition: RooAbsArg.h:609
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:302
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:544
void setAttribute(const Text_t *name, Bool_t value=kTRUE)
Set (default) or clear a named boolean attribute of this object.
Definition: RooAbsArg.cxx:257
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1718
std::map< std::string, std::string > _stringAttrib
Definition: RooAbsArg.h:610
RooAbsArg * findServer(const char *name) const
Return server of this with name name. Returns nullptr if not found.
Definition: RooAbsArg.h:203
OperMode operMode() const
Definition: RooAbsArg.h:474
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
Int_t getSize() const
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Storage_t::size_type size() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
const char * GetName() const
Returns name of object.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
friend class CacheElem
Definition: RooAbsPdf.h:333
virtual void resetErrorCounters(Int_t resetValue=10)
Reset error counter to given value, limiting the number of future error messages for this pdf to 'res...
Definition: RooAbsPdf.cxx:622
friend class RooRealIntegral
Definition: RooAbsPdf.h:313
Int_t _errorCount
Definition: RooAbsPdf.h:343
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:3283
RooArgSet * _normSet
Normalization integral (owned by _normMgr)
Definition: RooAbsPdf.h:321
static Int_t _verboseEval
Definition: RooAbsPdf.h:314
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Bool_t _forceNumInt
Definition: RooAbsReal.h:453
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, 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
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:560
Transiet cache with transformed values of coefficients.
Definition: RooAddModel.h:102
RooArgList _rangeProjList
Definition: RooAddModel.h:111
RooArgList _suppProjList
Definition: RooAddModel.h:109
RooArgList _refRangeProjList
Definition: RooAddModel.h:110
RooArgList _suppNormList
Definition: RooAddModel.h:104
virtual RooArgList containedArgs(Action)
List all RooAbsArg derived contents in this cache element.
virtual RooArgList containedArgs(Action)
List all RooAbsArg derived contents in this cache element.
Double_t evaluate() const
Calculate the current value.
friend class RooAddGenContext
Definition: RooAddModel.h:88
RooObjCacheManager _projCacheMgr
Definition: RooAddModel.h:116
Double_t analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const
Return analytical integral defined by given scenario code.
CacheElem * getProjCache(const RooArgSet *nset, const RooArgSet *iset=0, const char *rangeName=0) const
Retrieve cache element with for calculation of p.d.f value with normalization set nset and integrated...
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const
Return specialized context to efficiently generate toy events from RooAddPdfs.
void getCompIntList(const RooArgSet *nset, const RooArgSet *iset, pRooArgList &compIntList, Int_t &code, const char *isetRangeName) const
Check if this configuration was created before.
RooSetProxy _refCoefNorm
Definition: RooAddModel.h:95
virtual RooResolutionModel * convolution(RooFormulaVar *basis, RooAbsArg *owner) const
Instantiate a clone of this resolution model representing a convolution with given basis function.
virtual void selectNormalizationRange(const char *rangeName=0, Bool_t force=kFALSE)
Interface function used by test statistics to freeze choice of range for interpretation of fraction c...
Int_t _coefErrCount
Definition: RooAddModel.h:140
void updateCoefficients(CacheElem &cache, const RooArgSet *nset) const
Update the coefficient values in the given cache element: calculate new remainder fraction,...
RooArgSet _ownedComps
Coefficient error counter.
Definition: RooAddModel.h:142
virtual ~RooAddModel()
Destructor.
virtual void resetErrorCounters(Int_t resetValue=10)
Reset error counter to given value, limiting the number of future error messages for this pdf to 'res...
RooListProxy _coefList
Definition: RooAddModel.h:134
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=0) const
Variant of getAnalyticalIntegral that is also passed the normalization set that should be applied to ...
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return the number of expected events, which is either the sum of all coefficients or the sum of the c...
RooListProxy _pdfList
Registry of component analytical integration codes.
Definition: RooAddModel.h:133
RooObjCacheManager _intCacheMgr
Definition: RooAddModel.h:129
void fixCoefNormalization(const RooArgSet &refCoefNorm)
By default the interpretation of the fraction coefficients is performed in the contextual choice of o...
virtual void selectNormalization(const RooArgSet *depSet=0, Bool_t force=kFALSE)
Interface function used by test statistics to freeze choice of observables for interpretation of frac...
void generateEvent(Int_t code)
This function should never be called as RooAddModel implements a custom generator context.
Bool_t _projectCoefs
Reference range name for coefficient interpreation.
Definition: RooAddModel.h:98
void fixCoefRange(const char *rangeName)
By default the interpretation of the fraction coefficients is performed in the default range.
TNamed * _refCoefRangeName
Reference observable set for coefficient interpretation.
Definition: RooAddModel.h:96
virtual Bool_t checkObservables(const RooArgSet *nset) const
Check if PDF is valid for given normalization set.
void printMetaArgs(std::ostream &os) const
Customized printing of arguments of a RooAddModel to more intuitively reflect the contents of the pro...
Bool_t _haveLastCoef
List of supplemental normalization factors.
Definition: RooAddModel.h:137
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Return pseud-code that indicates if all components can do internal generation (1) or not (0)
Bool_t isDirectGenSafe(const RooAbsArg &arg) const
Direct generation is safe if all components say so.
Bool_t _allExtendable
Definition: RooAddModel.h:138
virtual Int_t basisCode(const char *name) const
Return code for basis function representing by 'name' string.
Double_t * _coefCache
Definition: RooAddModel.h:99
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
virtual Bool_t addOwned(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling addOwned() for each element in the source...
Definition: RooArgSet.h:92
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
const RooNameSet * nameSet2ByIndex(Int_t index) const
Retrieve RooNameSet associated with slot at given index.
T * getObjByIndex(Int_t index) const
Retrieve payload object by slot index.
void reset()
Clear the cache.
Int_t lastIndex() const
const RooNameSet * nameSet1ByIndex(Int_t index) const
Retrieve RooNameSet associated with slot at given index.
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition: RooFormulaVar.h:29
RooAbsArg * getParameter(const char *name) const
Definition: RooFormulaVar.h:40
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
static const char * str(const TNamed *ptr)
Return C++ string corresponding to given TNamed pointer.
Definition: RooNameReg.cxx:103
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:93
RooArgSet * select(const RooArgSet &list) const
Construct a RooArgSet of objects in input 'list' whose names match to those in the internal name list...
Definition: RooNameSet.cxx:185
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
static RooConstVar & value(Double_t value)
Return a constant value object with given value.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
RooResolutionModel is the base class for PDFs that represent a resolution model that can be convolute...
virtual void changeBasis(RooFormulaVar *basis)
Change the basis function we convolute with.
virtual Int_t basisCode(const char *name) const =0
virtual RooResolutionModel * convolution(RooFormulaVar *basis, RooAbsArg *owner) const
Instantiate a clone of this resolution model representing a convolution with given basis function.
RooFormulaVar * _basis
RooTemplateProxy< RooAbsRealLValue > x
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Overloaded RooArgSet::add() method inserts 'var' into set and registers 'var' as server to owner with...
virtual void removeAll()
Remove all argument inset using remove(const RooAbsArg&).
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
TString & Append(const char *cs)
Definition: TString.h:559
@ InputArguments
Definition: RooGlobalFunc.h:68
Definition: first.py:1