Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 _codeReg(10),
70 _snormList(0),
71 _haveLastCoef(false),
72 _allExtendable(false)
73{
74 _coefCache = new Double_t[10] ;
76}
77
78
79
80////////////////////////////////////////////////////////////////////////////////
81/// Generic constructor from list of PDFs and list of coefficients.
82/// Each pdf list element (i) is paired with coefficient list element (i).
83/// The number of coefficients must be either equal to the number of PDFs,
84/// in which case extended MLL fitting is enabled, or be one less.
85///
86/// All PDFs must inherit from RooAbsPdf. All coefficients must inherit from RooAbsReal.
87
88RooAddModel::RooAddModel(const char *name, const char *title, const RooArgList& inPdfList, const RooArgList& inCoefList, Bool_t ownPdfList) :
89 RooResolutionModel(name,title,((RooResolutionModel*)inPdfList.at(0))->convVar()),
90 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
91 _refCoefRangeName(0),
92 _projectCoefs(kFALSE),
93 _projCacheMgr(this,10),
94 _intCacheMgr(this,10),
95 _codeReg(10),
96 _pdfList("!pdfs","List of PDFs",this),
97 _coefList("!coefficients","List of coefficients",this),
98 _haveLastCoef(kFALSE),
99 _allExtendable(kFALSE)
100{
101 if (inPdfList.getSize()>inCoefList.getSize()+1) {
102 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName()
103 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
104 assert(0) ;
105 }
106
107 // Constructor with N PDFs and N or N-1 coefs
108 TIterator* pdfIter = inPdfList.createIterator() ;
109 TIterator* coefIter = inCoefList.createIterator() ;
110 RooAbsPdf* pdf ;
111 RooAbsReal* coef ;
112
113 while((coef = (RooAbsPdf*)coefIter->Next())) {
114 pdf = (RooAbsPdf*) pdfIter->Next() ;
115 if (!pdf) {
116 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName()
117 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
118 assert(0) ;
119 }
120 if (!dynamic_cast<RooAbsReal*>(coef)) {
121 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal, ignored" << endl ;
122 continue ;
123 }
124 if (!dynamic_cast<RooAbsReal*>(pdf)) {
125 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") pdf " << pdf->GetName() << " is not of type RooAbsPdf, ignored" << endl ;
126 continue ;
127 }
128 _pdfList.add(*pdf) ;
129 _coefList.add(*coef) ;
130 }
131
132 pdf = (RooAbsPdf*) pdfIter->Next() ;
133 if (pdf) {
134 if (!dynamic_cast<RooAbsReal*>(pdf)) {
135 coutE(InputArguments) << "RooAddModel::RooAddModel(" << GetName() << ") last pdf " << coef->GetName() << " is not of type RooAbsPdf, fatal error" << endl ;
136 assert(0) ;
137 }
138 _pdfList.add(*pdf) ;
139 } else {
141 }
142
143 delete pdfIter ;
144 delete coefIter ;
145
148
149 if (ownPdfList) {
151 }
152
153}
154
155
156
157////////////////////////////////////////////////////////////////////////////////
158/// Copy constructor
159
160RooAddModel::RooAddModel(const RooAddModel& other, const char* name) :
162 _refCoefNorm("!refCoefNorm",this,other._refCoefNorm),
163 _refCoefRangeName((TNamed*)other._refCoefRangeName),
164 _projectCoefs(other._projectCoefs),
165 _projCacheMgr(other._projCacheMgr,this),
166 _intCacheMgr(other._intCacheMgr,this),
167 _codeReg(other._codeReg),
168 _pdfList("!pdfs",this,other._pdfList),
169 _coefList("!coefficients",this,other._coefList),
170 _haveLastCoef(other._haveLastCoef),
171 _allExtendable(other._allExtendable)
172{
175}
176
177
178
179////////////////////////////////////////////////////////////////////////////////
180/// Destructor
181
183{
184 if (_coefCache) delete[] _coefCache ;
185}
186
187
188
189////////////////////////////////////////////////////////////////////////////////
190/// By default the interpretation of the fraction coefficients is
191/// performed in the contextual choice of observables. This makes the
192/// shape of the p.d.f explicitly dependent on the choice of
193/// observables. This method instructs RooAddModel to freeze the
194/// interpretation of the coefficients to be done in the given set of
195/// observables. If frozen, fractions are automatically transformed
196/// from the reference normalization set to the contextual normalization
197/// set by ratios of integrals
198
200{
201 if (refCoefNorm.getSize()==0) {
203 return ;
204 }
206
208 _refCoefNorm.add(refCoefNorm) ;
209
211}
212
213
214
215////////////////////////////////////////////////////////////////////////////////
216/// By default the interpretation of the fraction coefficients is
217/// performed in the default range. This make the shape of a RooAddModel
218/// explicitly dependent on the range of the observables. To allow
219/// a range independent definition of the fraction this function
220/// instructs RooAddModel to freeze its interpretation in the given
221/// named range. If the current normalization range is different
222/// from the reference range, the appropriate fraction coefficients
223/// are automically calculation from the reference fractions using
224/// ratios if integrals
225
226void RooAddModel::fixCoefRange(const char* rangeName)
227{
230}
231
232
233
234////////////////////////////////////////////////////////////////////////////////
235/// Instantiate a clone of this resolution model representing a convolution with given
236/// basis function. The owners object name is incorporated in the clones name
237/// to avoid multiple convolution objects with the same name in complex PDF structures.
238///
239/// RooAddModel will clone all the component models to create a composite convolution object
240
242{
243 // Check that primary variable of basis functions is our convolution variable
244 if (inBasis->getParameter(0) != x.absArg()) {
245 coutE(InputArguments) << "RooAddModel::convolution(" << GetName()
246 << ") convolution parameter of basis function and PDF don't match" << endl ;
247 ccoutE(InputArguments) << "basis->findServer(0) = " << inBasis->findServer(0) << " " << inBasis->findServer(0)->GetName() << endl ;
248 ccoutE(InputArguments) << "x.absArg() = " << x.absArg() << " " << x.absArg()->GetName() << endl ;
249 inBasis->Print("v") ;
250 return 0 ;
251 }
252
253 TString newName(GetName()) ;
254 newName.Append("_conv_") ;
255 newName.Append(inBasis->GetName()) ;
256 newName.Append("_[") ;
257 newName.Append(owner->GetName()) ;
258 newName.Append("]") ;
259
260 TString newTitle(GetTitle()) ;
261 newTitle.Append(" convoluted with basis function ") ;
262 newTitle.Append(inBasis->GetName()) ;
263
264 RooArgList modelList ;
265 for (auto obj : _pdfList) {
266 auto model = static_cast<RooResolutionModel*>(obj);
267 // Create component convolution
268 RooResolutionModel* conv = model->convolution(inBasis,owner) ;
269 modelList.add(*conv) ;
270 }
271
272 RooArgList theCoefList ;
273 for (auto coef : _coefList) {
274 theCoefList.add(*coef) ;
275 }
276
277 RooAddModel* convSum = new RooAddModel(newName,newTitle,modelList,theCoefList,kTRUE) ;
278 for (std::set<std::string>::const_iterator attrIt = _boolAttrib.begin();
279 attrIt != _boolAttrib.end(); ++attrIt) {
280 convSum->setAttribute((*attrIt).c_str()) ;
281 }
282 for (std::map<std::string,std::string>::const_iterator attrIt = _stringAttrib.begin();
283 attrIt != _stringAttrib.end(); ++attrIt) {
284 convSum->setStringAttribute((attrIt->first).c_str(), (attrIt->second).c_str()) ;
285 }
286 convSum->changeBasis(inBasis) ;
287 return convSum ;
288}
289
290
291
292////////////////////////////////////////////////////////////////////////////////
293/// Return code for basis function representing by 'name' string.
294/// The basis code of the first component model will be returned,
295/// if the basis is supported by all components. Otherwise 0
296/// is returned
297
299{
300 Bool_t first(kTRUE), code(0) ;
301 for (auto obj : _pdfList) {
302 auto model = static_cast<RooResolutionModel*>(obj);
303 Int_t subCode = model->basisCode(name) ;
304 if (first) {
305 code = subCode ;
306 first = kFALSE ;
307 } else if (subCode==0) {
308 code = 0 ;
309 }
310 }
311
312 return code ;
313}
314
315
316
317////////////////////////////////////////////////////////////////////////////////
318/// Retrieve cache element with for calculation of p.d.f value with normalization set nset and integrated over iset
319/// in range 'rangeName'. If cache element does not exist, create and fill it on the fly. The cache contains
320/// suplemental normalization terms (in case not all added p.d.f.s have the same observables), projection
321/// integrals to calculated transformed fraction coefficients when a frozen reference frame is provided
322/// and projection integrals for similar transformations when a frozen reference range is provided.
323
324RooAddModel::CacheElem* RooAddModel::getProjCache(const RooArgSet* nset, const RooArgSet* iset, const char* rangeName) const
325{
326 // Check if cache already exists
327 CacheElem* cache = (CacheElem*) _projCacheMgr.getObj(nset,iset,0,RooNameReg::ptr(rangeName)) ;
328 if (cache) {
329 return cache ;
330 }
331
332 //Create new cache
333 cache = new CacheElem ;
334
335 // *** PART 1 : Create supplemental normalization list ***
336
337 // Retrieve the combined set of dependents of this PDF ;
338 RooArgSet *fullDepList = getObservables(nset) ;
339 if (iset) {
340 fullDepList->remove(*iset,kTRUE,kTRUE) ;
341 }
342
343 // Fill with dummy unit RRVs for now
344 for (unsigned int i = 0; i < _pdfList.size(); ++i) {
345 auto pdf = static_cast<RooAbsPdf*>(&_pdfList[i]);
346 auto coef = i < _coefList.size() ? &_coefList[i] : nullptr;
347
348 // Start with full list of dependents
349 RooArgSet supNSet(*fullDepList) ;
350
351 // Remove PDF dependents
352 RooArgSet* pdfDeps = pdf->getObservables(nset) ;
353 if (pdfDeps) {
354 supNSet.remove(*pdfDeps,kTRUE,kTRUE) ;
355 delete pdfDeps ;
356 }
357
358 // Remove coef dependents
359 if (coef) {
360 RooArgSet* coefDeps = coef->getObservables(nset);
361 supNSet.remove(*coefDeps,kTRUE,kTRUE) ;
362 delete coefDeps ;
363 }
364
365 RooAbsReal* snorm ;
366 TString name(GetName()) ;
367 name.Append("_") ;
368 name.Append(pdf->GetName()) ;
369 name.Append("_SupNorm") ;
370 if (supNSet.getSize()>0) {
371 snorm = new RooRealIntegral(name,"Supplemental normalization integral",RooRealConstant::value(1.0),supNSet) ;
372 } else {
373 snorm = new RooRealVar(name,"Unit Supplemental normalization integral",1.0) ;
374 }
375 cache->_suppNormList.addOwned(*snorm) ;
376 }
377
378 delete fullDepList ;
379
380 if (_verboseEval>1) {
381 cxcoutD(Caching) << "RooAddModel::syncSuppNormList(" << GetName() << ") synching supplemental normalization list for norm" << (nset?*nset:RooArgSet()) << endl ;
382 if (dologD(Caching)) {
383 cache->_suppNormList.Print("v") ;
384 }
385 }
386
387
388 // *** PART 2 : Create projection coefficients ***
389
390 // If no projections required stop here
391 if (!_projectCoefs || _basis!=0 ) {
392 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
393 return cache ;
394 }
395
396
397 // Reduce iset/nset to actual dependents of this PDF
398 RooArgSet* nset2 = nset ? getObservables(nset) : new RooArgSet() ;
399
400 // Check if requested transformation is not identity
401 if (!nset2->equals(_refCoefNorm) || _refCoefRangeName !=0 || rangeName !=0 ) {
402
403 coutI(Caching) << "RooAddModel::syncCoefProjList(" << GetName() << ") creating coefficient projection integrals" << endl ;
404 ccoutI(Caching) << " from current normalization: " ; nset2->Print("1") ;
405 ccoutI(Caching) << " with current range: " << (rangeName?rangeName:"<none>") << endl ;
406 ccoutI(Caching) << " to reference normalization: " ; _refCoefNorm.Print("1") ;
407 ccoutI(Caching) << " with reference range: " << (_refCoefRangeName?RooNameReg::str(_refCoefRangeName):"<none>") << endl ;
408
409 // Recalculate projection integrals of PDFs
410 for (const auto obj : _pdfList) {
411 const auto thePdf = static_cast<RooAbsPdf*>(obj);
412
413 // Calculate projection integral
414 RooAbsReal* pdfProj ;
415 if (!nset2->equals(_refCoefNorm)) {
416 pdfProj = thePdf->createIntegral(*nset2,_refCoefNorm) ;
417 pdfProj->setOperMode(operMode()) ;
418 } else {
419 TString name(GetName()) ;
420 name.Append("_") ;
421 name.Append(thePdf->GetName()) ;
422 name.Append("_ProjectNorm") ;
423 pdfProj = new RooRealVar(name,"Unit Projection normalization integral",1.0) ;
424 }
425
426 cache->_projList.addOwned(*pdfProj) ;
427
428 // Calculation optional supplemental normalization term
429 RooArgSet supNormSet(_refCoefNorm) ;
430 RooArgSet* deps = thePdf->getParameters(RooArgSet()) ;
431 supNormSet.remove(*deps,kTRUE,kTRUE) ;
432 delete deps ;
433
434 RooAbsReal* snorm ;
435 TString name(GetName()) ;
436 name.Append("_") ;
437 name.Append(thePdf->GetName()) ;
438 name.Append("_ProjSupNorm") ;
439 if (supNormSet.getSize()>0) {
440 snorm = new RooRealIntegral(name,"Projection Supplemental normalization integral",
441 RooRealConstant::value(1.0),supNormSet) ;
442 } else {
443 snorm = new RooRealVar(name,"Unit Projection Supplemental normalization integral",1.0) ;
444 }
445 cache->_suppProjList.addOwned(*snorm) ;
446
447 // Calculate reference range adjusted projection integral
448 RooAbsReal* rangeProj1 ;
451 rangeProj1->setOperMode(operMode()) ;
452 } else {
453 TString theName(GetName()) ;
454 theName.Append("_") ;
455 theName.Append(thePdf->GetName()) ;
456 theName.Append("_RangeNorm1") ;
457 rangeProj1 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
458 }
459 cache->_refRangeProjList.addOwned(*rangeProj1) ;
460
461
462 // Calculate range adjusted projection integral
463 RooAbsReal* rangeProj2 ;
464 if (rangeName && _refCoefNorm.getSize()>0) {
465 rangeProj2 = thePdf->createIntegral(_refCoefNorm,_refCoefNorm,rangeName) ;
466 rangeProj2->setOperMode(operMode()) ;
467 } else {
468 TString theName(GetName()) ;
469 theName.Append("_") ;
470 theName.Append(thePdf->GetName()) ;
471 theName.Append("_RangeNorm2") ;
472 rangeProj2 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
473 }
474 cache->_rangeProjList.addOwned(*rangeProj2) ;
475
476 }
477
478 }
479
480 delete nset2 ;
481
482 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
483
484 return cache ;
485}
486
487
488
489////////////////////////////////////////////////////////////////////////////////
490/// Update the coefficient values in the given cache element: calculate new remainder
491/// fraction, normalize fractions obtained from extended ML terms to unity and
492/// multiply these the various range and dimensional corrections needed in the
493/// current use context
494
496{
497 // cxcoutD(ChangeTracking) << "RooAddModel::updateCoefficients(" << GetName() << ") update coefficients" << endl ;
498
499 Int_t i ;
500
501 // Straight coefficients
502 if (_allExtendable) {
503
504 // coef[i] = expectedEvents[i] / SUM(expectedEvents)
505 Double_t coefSum(0) ;
506 for (i=0 ; i<_pdfList.getSize() ; i++) {
508 coefSum += _coefCache[i] ;
509 }
510 if (coefSum==0.) {
511 coutW(Eval) << "RooAddModel::updateCoefCache(" << GetName() << ") WARNING: total number of expected events is 0" << endl ;
512 } else {
513 for (i=0 ; i<_pdfList.getSize() ; i++) {
514 _coefCache[i] /= coefSum ;
515 }
516 }
517
518 } else {
519 if (_haveLastCoef) {
520
521 // coef[i] = coef[i] / SUM(coef)
522 Double_t coefSum(0) ;
523 for (i=0 ; i<_coefList.getSize() ; i++) {
524 _coefCache[i] = ((RooAbsPdf*)_coefList.at(i))->getVal(nset) ;
525 coefSum += _coefCache[i] ;
526 }
527 for (i=0 ; i<_coefList.getSize() ; i++) {
528 _coefCache[i] /= coefSum ;
529 }
530 } else {
531
532 // coef[i] = coef[i] ; coef[n] = 1-SUM(coef[0...n-1])
533 Double_t lastCoef(1) ;
534 for (i=0 ; i<_coefList.getSize() ; i++) {
535 _coefCache[i] = ((RooAbsPdf*)_coefList.at(i))->getVal(nset) ;
536 cxcoutD(Caching) << "SYNC: orig coef[" << i << "] = " << _coefCache[i] << endl ;
537 lastCoef -= _coefCache[i] ;
538 }
539 _coefCache[_coefList.getSize()] = lastCoef ;
540 cxcoutD(Caching) << "SYNC: orig coef[" << _coefList.getSize() << "] = " << _coefCache[_coefList.getSize()] << endl ;
541
542
543 // Warn about coefficient degeneration
544 if ((lastCoef<-1e-05 || (lastCoef-1)>1e-5) && _coefErrCount-->0) {
545 coutW(Eval) << "RooAddModel::updateCoefCache(" << GetName()
546 << " WARNING: sum of PDF coefficients not in range [0-1], value="
547 << 1-lastCoef << endl ;
548 if (_coefErrCount==0) {
549 coutW(Eval) << " (no more will be printed)" << endl ;
550 }
551 }
552 }
553 }
554
555
556
557 // Stop here if not projection is required or needed
558 if ((!_projectCoefs) || cache._projList.getSize()==0) {
559 // cout << "SYNC no projection required rangeName = " << (rangeName?rangeName:"<none>") << endl ;
560 return ;
561 }
562
563 // Adjust coefficients for given projection
564 Double_t coefSum(0) ;
565 for (i=0 ; i<_pdfList.getSize() ; i++) {
566 GlobalSelectComponentRAII compRAII(true);
567
568 RooAbsReal* pp = ((RooAbsReal*)cache._projList.at(i)) ;
569 RooAbsReal* sn = ((RooAbsReal*)cache._suppProjList.at(i)) ;
570 RooAbsReal* r1 = ((RooAbsReal*)cache._refRangeProjList.at(i)) ;
571 RooAbsReal* r2 = ((RooAbsReal*)cache._rangeProjList.at(i)) ;
572
573 if (dologD(Eval)) {
574 cxcoutD(Eval) << "pp = " << pp->GetName() << endl
575 << "sn = " << sn->GetName() << endl
576 << "r1 = " << r1->GetName() << endl
577 << "r2 = " << r2->GetName() << endl ;
579 r1->printCompactTree(ccoutD(Eval)) ;
580 }
581
582 Double_t proj = pp->getVal()/sn->getVal()*(r2->getVal()/r1->getVal()) ;
583
584 _coefCache[i] *= proj ;
585 coefSum += _coefCache[i] ;
586 }
587 for (i=0 ; i<_pdfList.getSize() ; i++) {
588 _coefCache[i] /= coefSum ;
589// cout << "POST-SYNC coef[" << i << "] = " << _coefCache[i] << endl ;
590 }
591
592
593
594}
595
596
597
598////////////////////////////////////////////////////////////////////////////////
599/// Calculate the current value
600
602{
603 const RooArgSet* nset = _normSet ;
604 CacheElem* cache = getProjCache(nset) ;
605
606 updateCoefficients(*cache,nset) ;
607
608
609 // Do running sum of coef/pdf pairs, calculate lastCoef.
610 Double_t snormVal ;
611 Double_t value(0) ;
612 Int_t i(0) ;
613 for (auto obj : _pdfList) {
614 auto pdf = static_cast<RooAbsPdf*>(obj);
615
616 if (_coefCache[i]!=0.) {
617 snormVal = nset ? ((RooAbsReal*)cache->_suppNormList.at(i))->getVal() : 1.0 ;
618 Double_t pdfVal = pdf->getVal(nset) ;
619 // Double_t pdfNorm = pdf->getNorm(nset) ;
620 if (pdf->isSelectedComp()) {
621 value += pdfVal*_coefCache[i]/snormVal ;
622 cxcoutD(Eval) << "RooAddModel::evaluate(" << GetName() << ") value += ["
623 << pdf->GetName() << "] " << pdfVal << " * " << _coefCache[i] << " / " << snormVal << endl ;
624 }
625 }
626 i++ ;
627 }
628
629 return value ;
630}
631
632
633
634////////////////////////////////////////////////////////////////////////////////
635/// Reset error counter to given value, limiting the number
636/// of future error messages for this pdf to 'resetValue'
637
639{
641 _coefErrCount = resetValue ;
642}
643
644
645
646////////////////////////////////////////////////////////////////////////////////
647/// Check if PDF is valid for given normalization set.
648/// Coeffient and PDF must be non-overlapping, but pdf-coefficient
649/// pairs may overlap each other
650
652{
653 Bool_t ret(kFALSE) ;
654
655 for (unsigned int i = 0; i < _coefList.size(); ++i) {
656 auto pdf = &_pdfList[i];
657 auto coef = &_coefList[i];
658
659 if (pdf->observableOverlaps(nset,*coef)) {
660 coutE(InputArguments) << "RooAddModel::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName()
661 << " and PDF " << pdf->GetName() << " have one or more dependents in common" << endl ;
662 ret = kTRUE ;
663 }
664 }
665
666 return ret ;
667}
668
669
670
671////////////////////////////////////////////////////////////////////////////////
672
674 const RooArgSet* normSet, const char* rangeName) const
675{
676 if (_forceNumInt) return 0 ;
677
678 // Declare that we can analytically integrate all requested observables
679 analVars.add(allVars) ;
680
681 // Retrieve (or create) the required component integral list
682 Int_t code ;
683 RooArgList *cilist ;
684 getCompIntList(normSet,&allVars,cilist,code,rangeName) ;
685
686 return code+1 ;
687
688}
689
690
691
692////////////////////////////////////////////////////////////////////////////////
693/// Check if this configuration was created before
694
695void RooAddModel::getCompIntList(const RooArgSet* nset, const RooArgSet* iset, pRooArgList& compIntList, Int_t& code, const char* isetRangeName) const
696{
697 Int_t sterileIdx(-1) ;
698
699 IntCacheElem* cache = (IntCacheElem*) _intCacheMgr.getObj(nset,iset,&sterileIdx,RooNameReg::ptr(isetRangeName)) ;
700 if (cache) {
701 code = _intCacheMgr.lastIndex() ;
702 compIntList = &cache->_intList ;
703
704 return ;
705 }
706
707 // Create containers for partial integral components to be generated
708 cache = new IntCacheElem ;
709
710 // Fill Cache
711 for (auto obj : _pdfList) {
712 auto model = static_cast<RooResolutionModel*>(obj);
713
714 RooAbsReal* intPdf = model->createIntegral(*iset,nset,0,isetRangeName) ;
715 cache->_intList.addOwned(*intPdf) ;
716 }
717
718 // Store the partial integral list and return the assigned code ;
719 code = _intCacheMgr.setObj(nset,iset,(RooAbsCacheElement*)cache,RooNameReg::ptr(isetRangeName)) ;
720
721 // Fill references to be returned
722 compIntList = &cache->_intList ;
723}
724
725
726
727////////////////////////////////////////////////////////////////////////////////
728/// Return analytical integral defined by given scenario code
729
730Double_t RooAddModel::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
731{
732 // No integration scenario
733 if (code==0) {
734 return getVal(normSet) ;
735 }
736
737 // Partial integration scenarios
739
740 RooArgList* compIntList ;
741
742 // If cache has been sterilized, revive this slot
743 if (cache==0) {
745 RooArgSet* nset = _intCacheMgr.nameSet1ByIndex(code-1)->select(*vars) ;
746 RooArgSet* iset = _intCacheMgr.nameSet2ByIndex(code-1)->select(*vars) ;
747
748 Int_t code2(-1) ;
749 getCompIntList(nset,iset,compIntList,code2,rangeName) ;
750
751 delete vars ;
752 delete nset ;
753 delete iset ;
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)
#define ccoutE(a)
#define cxcoutD(a)
#define coutW(a)
#define dologD(a)
#define coutE(a)
#define ccoutI(a)
#define ccoutD(a)
const Bool_t kFALSE
Definition RtypesCore.h:92
const Bool_t kTRUE
Definition RtypesCore.h:91
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
@ kName
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:313
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.
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
std::set< std::string > _boolAttrib
Definition RooAbsArg.h:641
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:340
void setAttribute(const Text_t *name, Bool_t value=kTRUE)
Set (default) or clear a named boolean attribute of this object.
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Set the operation mode of this node.
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...
std::map< std::string, std::string > _stringAttrib
Definition RooAbsArg.h:642
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:502
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.
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...
RooArgSet const * _normSet
Normalization integral (owned by _normMgr)
Definition RooAbsPdf.h:322
friend class RooRealIntegral
Definition RooAbsPdf.h:314
Int_t _errorCount
Definition RooAbsPdf.h:352
static Int_t _verboseEval
Definition RooAbsPdf.h:315
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:478
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:91
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 ...
Transiet cache with transformed values of coefficients.
RooArgList _refRangeProjList
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:27
Double_t evaluate() const
Calculate the current value.
friend class RooAddGenContext
Definition RooAddModel.h:87
RooObjCacheManager _projCacheMgr
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 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:94
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
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.
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
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.
RooObjCacheManager _intCacheMgr
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:97
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:95
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.
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
virtual Int_t basisCode(const char *name) const
Return code for basis function representing by 'name' string.
Double_t * _coefCache
Definition RooAddModel.h:98
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:70
RooAbsArg * absArg() const
Definition RooArgProxy.h:37
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Add element to non-owning set.
Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE) override
Add element to an owning set.
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...
RooAbsArg * getParameter(const char *name) const
Return pointer to parameter with given name.
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.
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
RooArgSet * select(const RooArgSet &list) const
Construct a RooArgSet of objects in input 'list' whose names match to those in the internal name list...
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 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&).
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:136
TString & Append(const char *cs)
Definition TString.h:564
Definition first.py:1