Logo ROOT  
Reference Guide
RooAddPdf.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 RooAddPdf
19 \ingroup Roofitcore
20
21RooAddPdf is an efficient implementation of a sum of PDFs of the form
22
23\f[
24 \sum_{i=1}^{n} c_i \cdot \mathrm{PDF}_i
25\f]
26
27or
28\f[
29 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
30\f]
31
32The first form is for extended likelihood fits, where the
33expected number of events is \f$ \sum_i c_i \f$. The coefficients \f$ c_i \f$
34can either be explicitly provided, or, if all components support
35extended likelihood fits, they can be calculated from the contribution
36of each PDF to the total expected number of events.
37
38In the second form, the sum of the coefficients is required to be 1 or less,
39and the coefficient of the last PDF is calculated automatically from the condition
40that the sum of all coefficients has to be 1.
41
42### Recursive coefficients
43It is also possible to parameterise the coefficients recursively
44
45\f[
46 \sum_{i=1}^n c_i \prod_{j=1}^{i-1} \left[ (1-c_j) \right] \cdot \mathrm{PDF}_i \\
47 = c_1 \cdot \mathrm{PDF}_1 + (1-c_1)\, c_2 \cdot \mathrm{PDF}_2 + \ldots + (1-c_1)\ldots(1-c_{n-1}) \cdot 1 \cdot \mathrm{PDF}_n \\
48\f]
49
50In this form the sum of the coefficients is always less than 1.0
51for all possible values of the individual coefficients between 0 and 1.
52\note Don't pass the \f$ n^\mathrm{th} \f$ coefficient. It is always 1, since the normalisation condition removes one degree of freedom.
53
54RooAddPdf relies on each component PDF to be normalized and will perform
55no normalization other than calculating the proper last coefficient \f$ c_n \f$, if requested.
56An (enforced) condition for this assumption is that each \f$ \mathrm{PDF}_i \f$ is independent of each \f$ c_i \f$.
57
58## Difference between RooAddPdf / RooRealSumFunc / RooRealSumPdf
59- RooAddPdf is a PDF of PDFs, *i.e.* its components need to be normalised and non-negative.
60- RooRealSumPdf is a PDF of functions, *i.e.*, its components can be negative, but their sum cannot be. The normalisation
61 is computed automatically, unless the PDF is extended (see above).
62- RooRealSumFunc is a sum of functions. It is neither normalised, nor need it be positive.
63
64*/
65
66#include "RooAddPdf.h"
67
68#include "RooAddHelpers.h"
69#include "RooAddGenContext.h"
70#include "RooBatchCompute.h"
71#include "RooDataSet.h"
72#include "RooGlobalFunc.h"
73#include "RooNaNPacker.h"
74#include "RooRealProxy.h"
75#include "RooRealVar.h"
76#include "RooRealConstant.h"
77#include "RooRealIntegral.h"
78#include "RooRealSumPdf.h"
80
81#include <algorithm>
82#include <memory>
83#include <sstream>
84#include <set>
85
86using namespace std;
87
89
90
91////////////////////////////////////////////////////////////////////////////////
92/// Dummy constructor
93
94RooAddPdf::RooAddPdf(const char *name, const char *title) :
95 RooAbsPdf(name,title),
96 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,false,false),
97 _projCacheMgr(this,10),
98 _pdfList("!pdfs","List of PDFs",this),
99 _coefList("!coefficients","List of coefficients",this),
100 _coefErrCount{_errorCount}
101{
103}
104
105
107
108 // Two pdfs with the same name are only allowed in the input list if they are
109 // actually the same object.
110 using PdfInfo = std::pair<std::string,RooAbsArg*>;
111 std::set<PdfInfo> seen;
112 for(auto const& pdf : _pdfList) {
113 PdfInfo elem{pdf->GetName(), pdf};
114 auto comp = [&](PdfInfo const& p){ return p.first == elem.first && p.second != elem.second; };
115 auto found = std::find_if(seen.begin(), seen.end(), comp);
116 if(found != seen.end()) {
117 std::stringstream errorMsg;
118 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
119 << ") pdf list contains pdfs with duplicate name \"" << pdf->GetName() << "\"."
120 << std::endl;
121 coutE(InputArguments) << errorMsg.str();
122 throw std::invalid_argument(errorMsg.str().c_str());
123 }
124 seen.insert(elem);
125 }
126
127 _coefCache.resize(_pdfList.size());
128}
129
130
131////////////////////////////////////////////////////////////////////////////////
132/// Constructor with two PDFs and one coefficient
133
134RooAddPdf::RooAddPdf(const char *name, const char *title,
135 RooAbsPdf& pdf1, RooAbsPdf& pdf2, RooAbsReal& coef1) :
136 RooAddPdf(name, title)
137{
138 _pdfList.add(pdf1) ;
139 _pdfList.add(pdf2) ;
140 _coefList.add(coef1) ;
141
143}
144
145
146////////////////////////////////////////////////////////////////////////////////
147/// Generic constructor from list of PDFs and list of coefficients.
148/// Each pdf list element (i) is paired with coefficient list element (i).
149/// The number of coefficients must be either equal to the number of PDFs,
150/// in which case extended MLL fitting is enabled, or be one less.
151///
152/// All PDFs must inherit from RooAbsPdf. All coefficients must inherit from RooAbsReal
153///
154/// If the recursiveFraction flag is true, the coefficients are interpreted as recursive
155/// coefficients as explained in the class description.
156
157RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList, const RooArgList& inCoefList, bool recursiveFractions) :
158 RooAddPdf(name,title)
159{
160 _recursive = recursiveFractions;
161
162 if (inPdfList.size()>inCoefList.size()+1 || inPdfList.size()<inCoefList.size()) {
163 std::stringstream errorMsg;
164 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
165 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1." << endl ;
166 coutE(InputArguments) << errorMsg.str();
167 throw std::invalid_argument(errorMsg.str().c_str());
168 }
169
170 if (recursiveFractions && inPdfList.size()!=inCoefList.size()+1) {
171 std::stringstream errorMsg;
172 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
173 << "): Recursive fractions option can only be used if Npdf=Ncoef+1." << endl;
174 coutE(InputArguments) << errorMsg.str();
175 throw std::invalid_argument(errorMsg.str());
176 }
177
178 // Constructor with N PDFs and N or N-1 coefs
179 RooArgList partinCoefList ;
180
181 auto addRecursiveCoef = [this,&partinCoefList](RooAbsPdf& pdf, RooAbsReal& coef) -> RooAbsReal & {
182 partinCoefList.add(coef) ;
183 if(partinCoefList.size() == 1) {
184 // The first fraction is the first plain fraction
185 return coef;
186 }
187 // The i-th recursive fraction = (1-f1)*(1-f2)*...(fi) and is calculated from the list (f1,...,fi) by RooRecursiveFraction)
188 std::stringstream rfracName;
189 rfracName << GetName() << "_recursive_fraction_" << pdf.GetName() << "_" << partinCoefList.size();
190 auto rfrac = std::make_unique<RooRecursiveFraction>(rfracName.str().c_str(),"Recursive Fraction",partinCoefList) ;
191 auto & rfracRef = *rfrac;
192 addOwnedComponents(std::move(rfrac)) ;
193 return rfracRef;
194 };
195
196 for (auto i = 0u; i < inCoefList.size(); ++i) {
197 auto coef = dynamic_cast<RooAbsReal*>(inCoefList.at(i));
198 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(i));
199 if (inPdfList.at(i) == nullptr) {
200 std::stringstream errorMsg;
201 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
202 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
203 coutE(InputArguments) << errorMsg.str();
204 throw std::invalid_argument(errorMsg.str());
205 }
206 if (!coef) {
207 std::stringstream errorMsg;
208 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") coefficient " << (coef ? coef->GetName() : "") << " is not of type RooAbsReal, ignored" << endl ;
209 coutE(InputArguments) << errorMsg.str();
210 throw std::invalid_argument(errorMsg.str());
211 }
212 if (!pdf) {
213 std::stringstream errorMsg;
214 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "") << " is not of type RooAbsPdf, ignored" << endl ;
215 coutE(InputArguments) << errorMsg.str();
216 throw std::invalid_argument(errorMsg.str());
217 }
218 _pdfList.add(*pdf) ;
219
220 // Process recursive fraction mode separately
221 _coefList.add(recursiveFractions ? addRecursiveCoef(*pdf, *coef) : *coef);
222 }
223
224 if (inPdfList.size() == inCoefList.size() + 1) {
225 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(inCoefList.size()));
226
227 if (!pdf) {
228 coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") last argument " << inPdfList.at(inCoefList.size())->GetName() << " is not of type RooAbsPdf." << endl ;
229 throw std::invalid_argument("Last argument for RooAddPdf is not a PDF.");
230 }
231 _pdfList.add(*pdf) ;
232
233 // Process recursive fractions mode. Above, we verified that we don't have a last coefficient
234 if (recursiveFractions) {
235 _coefList.add(addRecursiveCoef(*pdf, RooFit::RooConst(1)));
236 // In recursive mode we always have Ncoef=Npdf, since we added it just above
237 _haveLastCoef=true ;
238 }
239
240 } else {
241 _haveLastCoef=true ;
242 }
243
245}
246
247
248////////////////////////////////////////////////////////////////////////////////
249/// Generic constructor from list of extended PDFs. There are no coefficients as the expected
250/// number of events from each components determine the relative weight of the PDFs.
251///
252/// All PDFs must inherit from RooAbsPdf.
253
254RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList) :
255 RooAddPdf(name,title)
256{
257 _allExtendable = true;
258
259 // Constructor with N PDFs
260 for (const auto pdfArg : inPdfList) {
261 auto pdf = dynamic_cast<const RooAbsPdf*>(pdfArg);
262
263 if (!pdf) {
264 coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "") << " is not of type RooAbsPdf, ignored" << endl ;
265 continue ;
266 }
267 if (!pdf->canBeExtended()) {
268 coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName() << " is not extendable, ignored" << endl ;
269 continue ;
270 }
271 _pdfList.add(*pdf) ;
272 }
273
275}
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// Copy constructor
280
281RooAddPdf::RooAddPdf(const RooAddPdf& other, const char* name) :
282 RooAbsPdf(other,name),
283 _refCoefNorm("!refCoefNorm",this,other._refCoefNorm),
284 _refCoefRangeName((TNamed*)other._refCoefRangeName),
285 _projectCoefs(other._projectCoefs),
286 _projCacheMgr(other._projCacheMgr,this),
287 _codeReg(other._codeReg),
288 _pdfList("!pdfs",this,other._pdfList),
289 _coefList("!coefficients",this,other._coefList),
290 _haveLastCoef(other._haveLastCoef),
291 _allExtendable(other._allExtendable),
292 _recursive(other._recursive)
293{
297}
298
299
300////////////////////////////////////////////////////////////////////////////////
301/// By default the interpretation of the fraction coefficients is
302/// performed in the contextual choice of observables. This makes the
303/// shape of the p.d.f explicitly dependent on the choice of
304/// observables. This method instructs RooAddPdf to freeze the
305/// interpretation of the coefficients to be done in the given set of
306/// observables. If frozen, fractions are automatically transformed
307/// from the reference normalization set to the contextual normalization
308/// set by ratios of integrals.
309
311{
312 if (refCoefNorm.empty()) {
313 _projectCoefs = false ;
314 return ;
315 }
316 _projectCoefs = true ;
317
319 _refCoefNorm.add(refCoefNorm) ;
320
322}
323
324
325
326////////////////////////////////////////////////////////////////////////////////
327/// By default, fraction coefficients are assumed to refer to the default
328/// fit range. This makes the shape of a RooAddPdf
329/// explicitly dependent on the range of the observables. Calling this function
330/// allows for a range-independent definition of the fractions, because it
331/// ties all coefficients to the given
332/// named range. If the normalisation range is different
333/// from this reference range, the appropriate fraction coefficients
334/// are automatically calculated from the reference fractions by
335/// integrating over the ranges, and comparing these integrals.
336
337void RooAddPdf::fixCoefRange(const char* rangeName)
338{
339 auto* newNamePtr = const_cast<TNamed*>(RooNameReg::ptr(rangeName));
340 if(newNamePtr != _refCoefRangeName) {
342 }
343 _refCoefRangeName = newNamePtr;
344 if (_refCoefRangeName) _projectCoefs = true ;
345}
346
347
348
349////////////////////////////////////////////////////////////////////////////////
350/// Retrieve cache element for the computation of the PDF normalisation.
351/// \param[in] nset Current normalisation set (integration over these variables yields 1).
352/// \param[in] iset Integration set. Variables to be integrated over (if integrations are performed).
353/// \param[in] rangeName Reference range for the integrals.
354///
355/// If a cache element does not exist, create and fill it on the fly. The cache also contains
356/// - Supplemental normalization terms (in case not all added p.d.f.s have the same observables)
357/// - Projection integrals to calculate transformed fraction coefficients when a frozen reference frame is provided
358/// - Projection integrals for similar transformations when a frozen reference range is provided.
359
360AddCacheElem* RooAddPdf::getProjCache(const RooArgSet* nset, const RooArgSet* iset, const char* rangeName) const
361{
362 // Check if cache already exists
363 auto cache = static_cast<AddCacheElem*>(_projCacheMgr.getObj(nset,iset,0,rangeName));
364 if (cache) {
365 return cache ;
366 }
367
368 //Create new cache
369 cache = new AddCacheElem{*this, _pdfList, _coefList, nset, iset, rangeName,
371
372 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
373
374 return cache;
375}
376
377
378////////////////////////////////////////////////////////////////////////////////
379/// Update the coefficient values in the given cache element: calculate new remainder
380/// fraction, normalize fractions obtained from extended ML terms to unity, and
381/// multiply the various range and dimensional corrections needed in the
382/// current use context.
383
385{
386 // Since this function updates the cache, it obviously needs write access:
387 auto& myCoefCache = const_cast<std::vector<double>&>(_coefCache);
388 myCoefCache.resize(_haveLastCoef ? _coefList.size() : _pdfList.size(), 0.);
389
390 // Straight coefficients
391 if (_allExtendable) {
392
393 // coef[i] = expectedEvents[i] / SUM(expectedEvents)
394 double coefSum(0) ;
395 std::size_t i = 0;
396 for (auto arg : _pdfList) {
397 auto pdf = static_cast<RooAbsPdf*>(arg);
398 myCoefCache[i] = pdf->expectedEvents(!_refCoefNorm.empty()?&_refCoefNorm:nset) ;
399 coefSum += myCoefCache[i] ;
400 i++ ;
401 }
402
403 if (coefSum==0.) {
404 coutW(Eval) << "RooAddPdf::updateCoefCache(" << GetName() << ") WARNING: total number of expected events is 0" << endl ;
405 } else {
406 for (std::size_t j=0; j < _pdfList.size(); j++) {
407 myCoefCache[j] /= coefSum ;
408 }
409 }
410
411 } else {
412 if (_haveLastCoef) {
413
414 // coef[i] = coef[i] / SUM(coef)
415 double coefSum(0) ;
416 std::size_t i=0;
417 for (auto coefArg : _coefList) {
418 auto coef = static_cast<RooAbsReal*>(coefArg);
419 myCoefCache[i] = coef->getVal(nset) ;
420 coefSum += myCoefCache[i++];
421 }
422 if (coefSum==0.) {
423 coutW(Eval) << "RooAddPdf::updateCoefCache(" << GetName() << ") WARNING: sum of coefficients is zero 0" << endl ;
424 } else {
425 const double invCoefSum = 1./coefSum;
426 for (std::size_t j=0; j < _coefList.size(); j++) {
427 myCoefCache[j] *= invCoefSum;
428 }
429 }
430 } else {
431
432 // coef[i] = coef[i] ; coef[n] = 1-SUM(coef[0...n-1])
433 double lastCoef(1) ;
434 std::size_t i=0;
435 for (auto coefArg : _coefList) {
436 auto coef = static_cast<RooAbsReal*>(coefArg);
437 myCoefCache[i] = coef->getVal(nset) ;
438 lastCoef -= myCoefCache[i++];
439 }
440 myCoefCache[_coefList.size()] = lastCoef ;
441
442 // Treat coefficient degeneration
443 const float coefDegen = lastCoef < 0. ? -lastCoef : (lastCoef > 1. ? lastCoef - 1. : 0.);
444 if (coefDegen > 1.E-5) {
445 myCoefCache[_coefList.size()] = RooNaNPacker::packFloatIntoNaN(100.f*coefDegen);
446
447 std::stringstream msg;
448 if (_coefErrCount-->0) {
449 msg << "RooAddPdf::updateCoefCache(" << GetName()
450 << " WARNING: sum of PDF coefficients not in range [0-1], value="
451 << 1-lastCoef ;
452 if (_coefErrCount==0) {
453 msg << " (no more will be printed)" ;
454 }
455 coutW(Eval) << msg.str() << std::endl;
456 }
457 }
458 }
459 }
460
461
462 // Stop here if not projection is required or needed
463 if ((!_projectCoefs && _normRange.Length()==0) || cache._projList.empty()) {
464 return ;
465 }
466
467 // Adjust coefficients for given projection
468 double coefSum(0) ;
469 {
471
472 for (std::size_t i = 0; i < _pdfList.size(); i++) {
473
474 RooAbsReal* pp = ((RooAbsReal*)cache._projList.at(i)) ;
475 RooAbsReal* sn = ((RooAbsReal*)cache._suppProjList.at(i)) ;
476 RooAbsReal* r1 = ((RooAbsReal*)cache._refRangeProjList.at(i)) ;
477 RooAbsReal* r2 = ((RooAbsReal*)cache._rangeProjList.at(i)) ;
478
479 double proj = pp->getVal()/sn->getVal()*(r2->getVal()/r1->getVal()) ;
480
481 myCoefCache[i] *= proj ;
482 coefSum += myCoefCache[i] ;
483 }
484 }
485
486
488 for (std::size_t i=0; i < _pdfList.size(); ++i) {
489 ccoutD(Caching) << " ALEX: POST-SYNC coef[" << i << "] = " << myCoefCache[i]
490 << " ( _coefCache[i]/coefSum = " << myCoefCache[i]*coefSum << "/" << coefSum << " ) "<< endl ;
491 }
492 }
493
494 if (coefSum==0.) {
495 coutE(Eval) << "RooAddPdf::updateCoefCache(" << GetName() << ") sum of coefficients is zero." << endl ;
496 }
497
498 for (std::size_t i=0; i < _pdfList.size(); i++) {
499 myCoefCache[i] /= coefSum ;
500 }
501
502}
503
504////////////////////////////////////////////////////////////////////////////////
505/// Look up projection cache and per-PDF norm sets. If a PDF doesn't have a special
506/// norm set, use the `defaultNorm`. If `defaultNorm == nullptr`, use the member
507/// _normSet.
508std::pair<const RooArgSet*, AddCacheElem*> RooAddPdf::getNormAndCache(const RooArgSet* nset) const {
509
510 // Treat empty normalization set and nullptr the same way.
511 if(nset && nset->empty()) nset = nullptr;
512
513 if (nset == nullptr) {
514 if (!_refCoefNorm.empty()) {
515 nset = &_refCoefNorm ;
516 }
517 }
518
519 // A RooAddPdf needs to have a normalization set defined, otherwise its
520 // coefficient will not be uniquely defined. Its shape depends on the
521 // normalization provided. Un-normalized calls to RooAddPdf can happen in
522 // Roofit, when printing the pdf's or when computing integrals. In these case,
523 // if the pdf has a normalization set previously defined (i.e. stored as a
524 // datamember in _copyOfLastNormSet) it should use it by default when the pdf
525 // is evaluated without passing a normalizations set (in pdf->getVal(nullptr) )
526 // In the case of no pre-defined normalization set exists, a warning will be
527 // produced, since the obtained value will be arbitrary. Note that to avoid
528 // unnecessary warning messages, when calling RooAbsPdf::printValue or
529 // RooAbsPdf::graphVizTree, the printing of the warning messages for the
530 // RooFit::Eval topic is explicitly disabled.
531 {
532 // If nset is still nullptr, get the pointer to a copy of the last-used
533 // normalization set. It nset is not nullptr, check whether the copy of
534 // the last-used normalization set needs an update.
535 if(nset == nullptr) {
536 nset = _copyOfLastNormSet.get();
537 } else if(nset->uniqueId() != _idOfLastUsedNormSet) {
538 _copyOfLastNormSet = std::make_unique<const RooArgSet>(*nset);
540 }
541
542 // If nset is STILL nullptr, print a warning.
543 if (nset == nullptr) {
544 coutW(Eval) << "Evaluating RooAddPdf without a defined normalization set. This can lead to ambiguos "
545 "coefficients definition and incorrect results."
546 << " Use RooAddPdf::fixCoefNormalization(nset) to provide a normalization set for "
547 "defining uniquely RooAddPdf coefficients!"
548 << std::endl;
549 }
550 }
551
552
553 AddCacheElem* cache = getProjCache(nset) ;
554 updateCoefficients(*cache,nset) ;
555
556 return {nset, cache};
557}
558
559
560////////////////////////////////////////////////////////////////////////////////
561/// Calculate and return the current value
562
563double RooAddPdf::getValV(const RooArgSet* normSet) const
564{
565 auto normAndCache = getNormAndCache(normSet);
566 const RooArgSet* nset = normAndCache.first;
567 AddCacheElem* cache = normAndCache.second;
568
569 // Process change in last data set used
570 bool nsetChanged(false) ;
572 nsetChanged = syncNormalization(nset) ;
573 }
574
575 // Do running sum of coef/pdf pairs, calculate lastCoef.
576 if (isValueDirty() || nsetChanged) {
577 _value = 0.0;
578
579 for (unsigned int i=0; i < _pdfList.size(); ++i) {
580 const auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
581 double snormVal = 1.;
582 if (cache->_needSupNorm) {
583 snormVal = ((RooAbsReal*)cache->_suppNormList.at(i))->getVal();
584 }
585
586 double pdfVal = pdf.getVal(nset);
587 if (pdf.isSelectedComp()) {
588 _value += pdfVal*_coefCache[i]/snormVal;
589 }
590 }
592 }
593
594 return _value;
595}
596
597
598////////////////////////////////////////////////////////////////////////////////
599/// Compute addition of PDFs in batches.
600void RooAddPdf::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
601{
604 AddCacheElem* cache = getNormAndCache(nullptr).second;
605 for (unsigned int pdfNo = 0; pdfNo < _pdfList.size(); ++pdfNo)
606 {
607 auto pdf = static_cast<RooAbsPdf*>(&_pdfList[pdfNo]);
608 if (pdf->isSelectedComp())
609 {
610 pdfs.push_back(dataMap.at(pdf));
611 coefs.push_back(_coefCache[pdfNo] / (cache->_needSupNorm ?
612 static_cast<RooAbsReal*>(cache->_suppNormList.at(pdfNo))->getVal() : 1) );
613 }
614 }
616 dispatch->compute(stream, RooBatchCompute::AddPdf, output, nEvents, pdfs, coefs);
617}
618
619
620////////////////////////////////////////////////////////////////////////////////
621/// Reset error counter to given value, limiting the number
622/// of future error messages for this pdf to 'resetValue'
623
625{
627 _coefErrCount = resetValue ;
628}
629
630
631
632////////////////////////////////////////////////////////////////////////////////
633/// Check if PDF is valid for given normalization set.
634/// Coeffient and PDF must be non-overlapping, but pdf-coefficient
635/// pairs may overlap each other
636
638{
640}
641
642
643////////////////////////////////////////////////////////////////////////////////
644/// Determine which part (if any) of given integral can be performed analytically.
645/// If any analytical integration is possible, return integration scenario code
646///
647/// RooAddPdf queries each component PDF for its analytical integration capability of the requested
648/// set ('allVars'). It finds the largest common set of variables that can be integrated
649/// by all components. If such a set exists, it reconfirms that each component is capable of
650/// analytically integrating the common set, and combines the components individual integration
651/// codes into a single integration code valid for RooAddPdf.
652
654 const RooArgSet* normSet, const char* rangeName) const
655{
656
657 RooArgSet allAnalVars(*std::unique_ptr<RooArgSet>{getObservables(allVars)}) ;
658
659 Int_t n(0) ;
660
661 // First iteration, determine what each component can integrate analytically
662 for (const auto pdfArg : _pdfList) {
663 auto pdf = static_cast<const RooAbsPdf *>(pdfArg);
664 RooArgSet subAnalVars ;
665 pdf->getAnalyticalIntegralWN(allVars,subAnalVars,normSet,rangeName) ;
666
667 // Observables that cannot be integrated analytically by this component are dropped from the common list
668 for (const auto arg : allVars) {
669 if (!subAnalVars.find(arg->GetName()) && pdf->dependsOn(*arg)) {
670 allAnalVars.remove(*arg,true,true) ;
671 }
672 }
673 n++ ;
674 }
675
676 // If no observables can be integrated analytically, return code 0 here
677 if (allAnalVars.empty()) {
678 return 0 ;
679 }
680
681
682 // Now retrieve codes for integration over common set of analytically integrable observables for each component
683 n=0 ;
684 std::vector<Int_t> subCode(_pdfList.size());
685 bool allOK(true) ;
686 for (const auto arg : _pdfList) {
687 auto pdf = static_cast<const RooAbsPdf *>(arg);
688 RooArgSet subAnalVars ;
689 auto allAnalVars2 = std::unique_ptr<RooArgSet>{pdf->getObservables(allAnalVars)} ;
690 subCode[n] = pdf->getAnalyticalIntegralWN(*allAnalVars2,subAnalVars,normSet,rangeName) ;
691 if (subCode[n]==0 && !allAnalVars2->empty()) {
692 coutE(InputArguments) << "RooAddPdf::getAnalyticalIntegral(" << GetName() << ") WARNING: component PDF " << pdf->GetName()
693 << " advertises inconsistent set of integrals (e.g. (X,Y) but not X or Y individually."
694 << " Distributed analytical integration disabled. Please fix PDF" << endl ;
695 allOK = false ;
696 }
697 n++ ;
698 }
699 if (!allOK) {
700 return 0 ;
701 }
702
703 // Mare all analytically integrated observables as such
704 analVars.add(allAnalVars) ;
705
706 // Store set of variables analytically integrated
707 RooArgSet* intSet = new RooArgSet(allAnalVars) ;
708 Int_t masterCode = _codeReg.store(subCode,intSet)+1 ;
709
710 return masterCode ;
711}
712
713
714
715////////////////////////////////////////////////////////////////////////////////
716/// Return analytical integral defined by given scenario code
717
718double RooAddPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
719{
720 // WVE needs adaptation to handle new rangeName feature
721 if (code==0) {
722 return getVal(normSet) ;
723 }
724
725 // Retrieve analytical integration subCodes and set of observabels integrated over
726 RooArgSet* intSet ;
727 const std::vector<Int_t>& subCode = _codeReg.retrieve(code-1,intSet) ;
728 if (subCode.empty()) {
729 coutE(InputArguments) << "RooAddPdf::analyticalIntegral(" << GetName() << "): ERROR unrecognized integration code, " << code << endl ;
730 assert(0) ;
731 }
732
733 cxcoutD(Caching) << "RooAddPdf::aiWN(" << GetName() << ") calling getProjCache with nset = " << (normSet?*normSet:RooArgSet()) << endl ;
734
735 if ((normSet==0 || normSet->empty()) && !_refCoefNorm.empty()) {
736// cout << "WVE integration of RooAddPdf without normalization, but have reference set, using ref set for normalization" << endl ;
737 normSet = &_refCoefNorm ;
738 }
739
740 AddCacheElem* cache = getProjCache(normSet,intSet,0) ; // WVE rangename here?
741 updateCoefficients(*cache,normSet) ;
742
743 // Calculate the current value of this object
744 double value(0) ;
745
746 // Do running sum of coef/pdf pairs, calculate lastCoef.
747 double snormVal ;
748
749 //cout << "ROP::aIWN updateCoefCache with rangeName = " << (rangeName?rangeName:"<null>") << endl ;
750 RooArgList* snormSet = (!cache->_suppNormList.empty()) ? &cache->_suppNormList : nullptr ;
751 for (std::size_t i = 0; i < _pdfList.size(); ++i ) {
752 auto pdf = static_cast<const RooAbsPdf*>(_pdfList.at(i));
753
754 if (_coefCache[i]) {
755 snormVal = snormSet ? ((RooAbsReal*) cache->_suppNormList.at(i))->getVal() : 1.0 ;
756
757 // WVE swap this?
758 double val = pdf->analyticalIntegralWN(subCode[i],normSet,rangeName) ;
759 if (pdf->isSelectedComp()) {
760 value += val*_coefCache[i]/snormVal ;
761 }
762 }
763 }
764
765 return value ;
766}
767
768
769
770////////////////////////////////////////////////////////////////////////////////
771/// Return the number of expected events, which is either the sum of all coefficients
772/// or the sum of the components extended terms, multiplied with the fraction that
773/// is in the current range w.r.t the reference range
774
775double RooAddPdf::expectedEvents(const RooArgSet* nset) const
776{
777 double expectedTotal{0.0};
778
779 cxcoutD(Caching) << "RooAddPdf::expectedEvents(" << GetName() << ") calling getProjCache with nset = " << (nset?*nset:RooArgSet()) << endl ;
780 AddCacheElem& cache = *getProjCache(nset) ;
781 updateCoefficients(cache,nset) ;
782
783 if (!cache._rangeProjList.empty()) {
784
785 for (std::size_t i = 0; i < _pdfList.size(); ++i) {
786 auto const& r1 = static_cast<RooAbsReal&>(cache._refRangeProjList[i]);
787 auto const& r2 = static_cast<RooAbsReal&>(cache._rangeProjList[i]);
788 double ncomp = _allExtendable ? static_cast<RooAbsPdf&>(_pdfList[i]).expectedEvents(nset)
789 : static_cast<RooAbsReal&>(_coefList[i]).getVal(nset);
790 expectedTotal += (r2.getVal()/r1.getVal()) * ncomp ;
791
792 }
793
794 } else {
795
796 if (_allExtendable) {
797 for(auto const& arg : _pdfList) {
798 expectedTotal += static_cast<RooAbsPdf*>(arg)->expectedEvents(nset) ;
799 }
800 } else {
801 for(auto const& arg : _coefList) {
802 expectedTotal += static_cast<RooAbsReal*>(arg)->getVal(nset) ;
803 }
804 }
805
806 }
807 return expectedTotal ;
808}
809
810
811
812////////////////////////////////////////////////////////////////////////////////
813/// Interface function used by test statistics to freeze choice of observables
814/// for interpretation of fraction coefficients
815
816void RooAddPdf::selectNormalization(const RooArgSet* depSet, bool force)
817{
818
819 if (!force && !_refCoefNorm.empty()) {
820 return ;
821 }
822
823 if (!depSet) {
825 return ;
826 }
827
828 fixCoefNormalization(*std::unique_ptr<RooArgSet>{getObservables(depSet)}) ;
829}
830
831
832
833////////////////////////////////////////////////////////////////////////////////
834/// Interface function used by test statistics to freeze choice of range
835/// for interpretation of fraction coefficients
836
837void RooAddPdf::selectNormalizationRange(const char* rangeName, bool force)
838{
839 if (!force && _refCoefRangeName) {
840 return ;
841 }
842
843 fixCoefRange(rangeName) ;
844}
845
846
847
848////////////////////////////////////////////////////////////////////////////////
849/// Return specialized context to efficiently generate toy events from RooAddPdfs
850/// return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ; // WVE DEBUG
851
853 const RooArgSet* auxProto, bool verbose) const
854{
855 return RooAddGenContext::create(*this,vars,prototype,auxProto,verbose).release();
856}
857
858
859
860////////////////////////////////////////////////////////////////////////////////
861/// Loop over components for plot sampling hints and merge them if there are multiple
862
863std::list<double>* RooAddPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
864{
865 return RooRealSumPdf::plotSamplingHint(_pdfList, obs, xlo, xhi);
866}
867
868
869////////////////////////////////////////////////////////////////////////////////
870/// Loop over components for plot sampling hints and merge them if there are multiple
871
872std::list<double>* RooAddPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
873{
874 return RooRealSumPdf::binBoundaries(_pdfList, obs, xlo, xhi);
875}
876
877
878////////////////////////////////////////////////////////////////////////////////
879/// If all components that depend on obs are binned, so is their sum.
881{
883}
884
885
886////////////////////////////////////////////////////////////////////////////////
887/// Label OK'ed components of a RooAddPdf with cache-and-track
888
890{
892}
893
894
895
896////////////////////////////////////////////////////////////////////////////////
897/// Customized printing of arguments of a RooAddPdf to more intuitively reflect the contents of the
898/// product operator construction
899
900void RooAddPdf::printMetaArgs(std::ostream& os) const
901{
903}
#define cxcoutD(a)
Definition: RooMsgService.h:85
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define ccoutD(a)
Definition: RooMsgService.h:41
#define TRACE_CREATE
Definition: RooTrace.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
winID h TVirtualViewer3D TVirtualGLPainter p
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
RooArgList _projList
Projection integrals to be multiplied with coefficients.
Definition: RooAddHelpers.h:31
RooArgList _suppProjList
Projection integrals to multiply with coefficients for supplemental normalization.
Definition: RooAddHelpers.h:32
RooArgList _refRangeProjList
Range integrals to be multiplied with coefficients (reference range)
Definition: RooAddHelpers.h:33
RooArgList _suppNormList
Supplemental normalization list.
Definition: RooAddHelpers.h:28
RooArgList _rangeProjList
Range integrals to be multiplied with coefficients (target range)
Definition: RooAddHelpers.h:34
bool _needSupNorm
Does the above list contain any non-unit entries?
Definition: RooAddHelpers.h:29
const std::vector< Int_t > & retrieve(Int_t masterCode) const
Retrieve the array of integer codes associated with the given master code.
Int_t store(const std::vector< Int_t > &codeList, RooArgSet *set1=0, RooArgSet *set2=0, RooArgSet *set3=0, RooArgSet *set4=0)
Store given arrays of integer codes, and up to four RooArgSets in the registry (each setX pointer may...
void clearValueAndShapeDirty() const
Definition: RooAbsArg.h:611
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:805
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:308
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2163
bool isValueDirty() const
Definition: RooAbsArg.h:436
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool empty() const
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
virtual bool syncNormalization(const RooArgSet *dset, bool adjustProxies=true) const
Verify that the normalization integral cached with this PDF is valid for given set of normalization o...
Definition: RooAbsPdf.cxx:538
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:631
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition: RooAbsPdf.cxx:3195
TString _normRange
Normalization range.
Definition: RooAbsPdf.h:378
RooArgSet const * _normSet
Normalization integral (owned by _normMgr)
Definition: RooAbsPdf.h:343
RooAbsReal * _norm
Definition: RooAbsPdf.h:342
Int_t _errorCount
Number of errors remaining to print.
Definition: RooAbsPdf.h:370
static Int_t _verboseEval
Definition: RooAbsPdf.h:337
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
double _value
Cache for current value of object.
Definition: RooAbsReal.h:480
static std::unique_ptr< RooAbsGenContext > create(const Pdf_t &pdf, const RooArgSet &vars, const RooDataSet *prototype, const RooArgSet *auxProto, bool verbose)
Returns a RooAddGenContext if possible, or, if the RooAddGenContext doesn't support this particular R...
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
Definition: RooAddPdf.h:34
RooListProxy _coefList
List of coefficients.
Definition: RooAddPdf.h:129
bool _allExtendable
Flag indicating if all PDF components are extendable.
Definition: RooAddPdf.h:133
RooAICRegistry _codeReg
! Registry of component analytical integration codes
Definition: RooAddPdf.h:126
RooFit::UniqueId< RooArgSet >::Value_t _idOfLastUsedNormSet
!
Definition: RooAddPdf.h:147
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Return analytical integral defined by given scenario code.
Definition: RooAddPdf.cxx:718
std::unique_ptr< const RooArgSet > _copyOfLastNormSet
!
Definition: RooAddPdf.h:148
Int_t _coefErrCount
! Coefficient error counter
Definition: RooAddPdf.h:136
bool _haveLastCoef
Flag indicating if last PDFs coefficient was supplied in the ctor.
Definition: RooAddPdf.h:132
void selectNormalization(const RooArgSet *depSet=nullptr, bool force=false) override
Interface function used by test statistics to freeze choice of observables for interpretation of frac...
Definition: RooAddPdf.cxx:816
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooAddPdf to more intuitively reflect the contents of the produ...
Definition: RooAddPdf.cxx:900
void finalizeConstruction()
Definition: RooAddPdf.cxx:106
void setCacheAndTrackHints(RooArgSet &) override
Label OK'ed components of a RooAddPdf with cache-and-track.
Definition: RooAddPdf.cxx:889
bool _recursive
Flag indicating is fractions are treated recursively.
Definition: RooAddPdf.h:134
RooObjCacheManager _projCacheMgr
Definition: RooAddPdf.h:107
RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const override
Return specialized context to efficiently generate toy events from RooAddPdfs return RooAbsPdf::genCo...
Definition: RooAddPdf.cxx:852
bool checkObservables(const RooArgSet *nset) const override
Check if PDF is valid for given normalization set.
Definition: RooAddPdf.cxx:637
RooAddPdf()
Definition: RooAddPdf.h:37
void fixCoefNormalization(const RooArgSet &refCoefNorm)
By default the interpretation of the fraction coefficients is performed in the contextual choice of o...
Definition: RooAddPdf.cxx:310
std::pair< const RooArgSet *, AddCacheElem * > getNormAndCache(const RooArgSet *nset) const
Look up projection cache and per-PDF norm sets.
Definition: RooAddPdf.cxx:508
void updateCoefficients(AddCacheElem &cache, const RooArgSet *nset) const
Update the coefficient values in the given cache element: calculate new remainder fraction,...
Definition: RooAddPdf.cxx:384
RooSetProxy _refCoefNorm
Reference observable set for coefficient interpretation.
Definition: RooAddPdf.h:100
void selectNormalizationRange(const char *rangeName=nullptr, bool force=false) override
Interface function used by test statistics to freeze choice of range for interpretation of fraction c...
Definition: RooAddPdf.cxx:837
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Determine which part (if any) of given integral can be performed analytically.
Definition: RooAddPdf.cxx:653
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...
Definition: RooAddPdf.cxx:624
double expectedEvents(const RooArgSet *nset) const override
Return expected number of events for extended likelihood calculation, which is the sum of all coeffic...
Definition: RooAddPdf.cxx:775
double getValV(const RooArgSet *set=nullptr) const override
Calculate and return the current value.
Definition: RooAddPdf.cxx:563
void fixCoefRange(const char *rangeName)
By default, fraction coefficients are assumed to refer to the default fit range.
Definition: RooAddPdf.cxx:337
std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override
Loop over components for plot sampling hints and merge them if there are multiple.
Definition: RooAddPdf.cxx:863
RooListProxy _pdfList
List of component PDFs.
Definition: RooAddPdf.h:128
TNamed * _refCoefRangeName
Reference range name for coefficient interpreation.
Definition: RooAddPdf.h:101
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute addition of PDFs in batches.
Definition: RooAddPdf.cxx:600
bool isBinnedDistribution(const RooArgSet &obs) const override
If all components that depend on obs are binned, so is their sum.
Definition: RooAddPdf.cxx:880
AddCacheElem * getProjCache(const RooArgSet *nset, const RooArgSet *iset=nullptr, const char *rangeName=nullptr) const
Manager of cache with coefficient projections and transformations.
Definition: RooAddPdf.cxx:360
bool _projectCoefs
If true coefficients need to be projected for use in evaluate()
Definition: RooAddPdf.h:103
std::vector< double > _coefCache
! Transient cache with transformed values of coefficients
Definition: RooAddPdf.h:104
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Loop over components for plot sampling hints and merge them if there are multiple.
Definition: RooAddPdf.cxx:872
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:56
RooFit::UniqueId< RooArgSet > const & uniqueId() const
Returns a unique ID that is different for every instantiated RooArgSet.
Definition: RooArgSet.h:192
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
void reset()
Clear the cache.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter 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
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
static RooMsgService & instance()
Return reference to singleton instance.
static Int_t _debugCount
bool isActive(const RooAbsArg *self, RooFit::MsgTopic facility, RooFit::MsgLevel level)
Check if logging is active for given object/topic/RooFit::MsgLevel combination.
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:81
void setCacheAndTrackHints(RooArgSet &) override
Label OK'ed components of a RooRealSumPdf with cache-and-track.
bool checkObservables(const RooArgSet *nset) const override
Check if FUNC is valid for given normalization set.
std::list< double > * plotSamplingHint(RooAbsRealLValue &, double, double) const override
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Retrieve bin boundaries if this distribution is binned in obs.
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooRealSumPdf to more intuitively reflect the contents of the p...
bool isBinnedDistribution(const RooArgSet &obs) const override
Check if all components that depend on obs are binned.
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
Ssiz_t Length() const
Definition: TString.h:410
RooConstVar & RooConst(double val)
const Int_t n
Definition: legend1.C:16
for(Int_t i=0;i< n;i++)
Definition: legend1.C:18
std::vector< RooSpan< const double > > VarVector
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
std::vector< double > ArgVector
@ InputArguments
Definition: RooGlobalFunc.h:64
UniqueId< Class > const & getUniqueId(Class const *ptr)
A helper function to replace pointer comparisons with UniqueId comparisons.
Definition: UniqueId.h:87
static double packFloatIntoNaN(float payload)
Pack float into mantissa of a NaN.
Definition: RooNaNPacker.h:109
static void output()