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