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