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 "RooRealProxy.h"
74#include "RooRealVar.h"
75#include "RooRealConstant.h"
76#include "RooRealSumPdf.h"
78
79#include <algorithm>
80#include <memory>
81#include <sstream>
82#include <set>
83
85
86
87////////////////////////////////////////////////////////////////////////////////
88/// Dummy constructor
89
90RooAddPdf::RooAddPdf(const char *name, const char *title) :
91 RooAbsPdf(name,title),
92 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,false,false),
93 _projCacheMgr(this,10),
94 _pdfList("!pdfs","List of PDFs",this),
95 _coefList("!coefficients","List of coefficients",this),
96 _coefErrCount{_errorCount}
97{
99}
100
101
103
104 // Two pdfs with the same name are only allowed in the input list if they are
105 // actually the same object.
106 using PdfInfo = std::pair<std::string,RooAbsArg*>;
107 std::set<PdfInfo> seen;
108 for(auto const& pdf : _pdfList) {
109 PdfInfo elem{pdf->GetName(), pdf};
110 auto comp = [&](PdfInfo const& p){ return p.first == elem.first && p.second != elem.second; };
111 auto found = std::find_if(seen.begin(), seen.end(), comp);
112 if(found != seen.end()) {
113 std::stringstream errorMsg;
114 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
115 << ") pdf list contains pdfs with duplicate name \"" << pdf->GetName() << "\".";
116 coutE(InputArguments) << errorMsg.str() << std::endl;
117 throw std::invalid_argument(errorMsg.str().c_str());
118 }
119 seen.insert(elem);
120 }
121}
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Constructor with two PDFs and one coefficient
126
127RooAddPdf::RooAddPdf(const char *name, const char *title,
128 RooAbsPdf& pdf1, RooAbsPdf& pdf2, RooAbsReal& coef1) :
129 RooAddPdf(name, title)
130{
131 _pdfList.add(pdf1) ;
132 _pdfList.add(pdf2) ;
133 _coefList.add(coef1) ;
134
136}
137
138
139////////////////////////////////////////////////////////////////////////////////
140/// Generic constructor from list of PDFs and list of coefficients.
141/// Each pdf list element (i) is paired with coefficient list element (i).
142/// The number of coefficients must be either equal to the number of PDFs,
143/// in which case extended MLL fitting is enabled, or be one less.
144///
145/// All PDFs must inherit from RooAbsPdf. All coefficients must inherit from RooAbsReal
146///
147/// If the recursiveFraction flag is true, the coefficients are interpreted as recursive
148/// coefficients as explained in the class description.
149
150RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList, const RooArgList& inCoefList, bool recursiveFractions) :
151 RooAddPdf(name,title)
152{
153 _recursive = recursiveFractions;
154
155 if (inPdfList.size()>inCoefList.size()+1 || inPdfList.size()<inCoefList.size()) {
156 std::stringstream errorMsg;
157 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
158 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1.";
159 coutE(InputArguments) << errorMsg.str() << std::endl;
160 throw std::invalid_argument(errorMsg.str().c_str());
161 }
162
163 if (recursiveFractions && inPdfList.size()!=inCoefList.size()+1) {
164 std::stringstream errorMsg;
165 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
166 << "): Recursive fractions option can only be used if Npdf=Ncoef+1.";
167 coutE(InputArguments) << errorMsg.str() << std::endl;
168 throw std::invalid_argument(errorMsg.str());
169 }
170
171 // Constructor with N PDFs and N or N-1 coefs
172 RooArgList partinCoefList ;
173
174 auto addRecursiveCoef = [this,&partinCoefList](RooAbsPdf& pdf, RooAbsReal& coef) -> RooAbsReal & {
175 partinCoefList.add(coef) ;
176 if(partinCoefList.size() == 1) {
177 // The first fraction is the first plain fraction
178 return coef;
179 }
180 // The i-th recursive fraction = (1-f1)*(1-f2)*...(fi) and is calculated from the list (f1,...,fi) by RooRecursiveFraction)
181 std::stringstream rfracName;
182 rfracName << GetName() << "_recursive_fraction_" << pdf.GetName() << "_" << partinCoefList.size();
183 auto rfrac = std::make_unique<RooRecursiveFraction>(rfracName.str().c_str(),"Recursive Fraction",partinCoefList) ;
184 auto & rfracRef = *rfrac;
185 addOwnedComponents(std::move(rfrac)) ;
186 return rfracRef;
187 };
188
189 for (auto i = 0u; i < inCoefList.size(); ++i) {
190 auto coef = dynamic_cast<RooAbsReal*>(inCoefList.at(i));
191 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(i));
192 if (inPdfList.at(i) == nullptr) {
193 std::stringstream errorMsg;
194 errorMsg << "RooAddPdf::RooAddPdf(" << GetName()
195 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1";
196 coutE(InputArguments) << errorMsg.str() << std::endl;
197 throw std::invalid_argument(errorMsg.str());
198 }
199 if (!coef) {
200 std::stringstream errorMsg;
201 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") coefficient " << (coef ? coef->GetName() : "") << " is not of type RooAbsReal, ignored";
202 coutE(InputArguments) << errorMsg.str() << std::endl;
203 throw std::invalid_argument(errorMsg.str());
204 }
205 if (!pdf) {
206 std::stringstream errorMsg;
207 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "") << " is not of type RooAbsPdf, ignored";
208 coutE(InputArguments) << errorMsg.str() << std::endl;
209 throw std::invalid_argument(errorMsg.str());
210 }
211 _pdfList.add(*pdf) ;
212
213 // Process recursive fraction mode separately
214 _coefList.add(recursiveFractions ? addRecursiveCoef(*pdf, *coef) : *coef);
215 }
216
217 if (inPdfList.size() == inCoefList.size() + 1) {
218 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(inCoefList.size()));
219
220 if (!pdf) {
221 coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") last argument " << inPdfList.at(inCoefList.size())->GetName() << " is not of type RooAbsPdf." << std::endl;
222 throw std::invalid_argument("Last argument for RooAddPdf is not a PDF.");
223 }
224 _pdfList.add(*pdf) ;
225
226 // Process recursive fractions mode. Above, we verified that we don't have a last coefficient
227 if (recursiveFractions) {
228 _coefList.add(addRecursiveCoef(*pdf, RooFit::RooConst(1)));
229 // In recursive mode we always have Ncoef=Npdf, since we added it just above
230 _haveLastCoef=true ;
231 }
232
233 } else {
234 _haveLastCoef=true ;
235 }
236
238}
239
240
241////////////////////////////////////////////////////////////////////////////////
242/// Generic constructor from list of extended PDFs. There are no coefficients as the expected
243/// number of events from each components determine the relative weight of the PDFs.
244///
245/// All PDFs must inherit from RooAbsPdf.
246
247RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList) :
248 RooAddPdf(name,title)
249{
250 _allExtendable = true;
251
252 // Constructor with N PDFs
253 for (const auto pdfArg : inPdfList) {
254 auto pdf = dynamic_cast<const RooAbsPdf*>(pdfArg);
255
256 if (!pdf) {
257 std::stringstream errorMsg;
258 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "")
259 << " is not of type RooAbsPdf, RooAddPdf constructor call is invalid!";
260 coutE(InputArguments) << errorMsg.str() << std::endl;
261 throw std::invalid_argument(errorMsg.str().c_str());
262 }
263 if (!pdf->canBeExtended()) {
264 std::stringstream errorMsg;
265 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName()
266 << " is not extendable, RooAddPdf constructor call is invalid!";
267 coutE(InputArguments) << errorMsg.str() << std::endl;
268 throw std::invalid_argument(errorMsg.str().c_str());
269 }
270 _pdfList.add(*pdf) ;
271 }
272
274}
275
276
277////////////////////////////////////////////////////////////////////////////////
278/// Copy constructor
279
280RooAddPdf::RooAddPdf(const RooAddPdf& other, const char* name) :
281 RooAbsPdf(other,name),
282 _refCoefNorm("!refCoefNorm",this,other._refCoefNorm),
283 _refCoefRangeName((TNamed*)other._refCoefRangeName),
284 _projCacheMgr(other._projCacheMgr,this),
285 _codeReg(other._codeReg),
286 _pdfList("!pdfs",this,other._pdfList),
287 _coefList("!coefficients",this,other._coefList),
288 _haveLastCoef(other._haveLastCoef),
289 _allExtendable(other._allExtendable),
290 _recursive(other._recursive)
291{
295}
296
297
298////////////////////////////////////////////////////////////////////////////////
299/// By default the interpretation of the fraction coefficients is
300/// performed in the contextual choice of observables. This makes the
301/// shape of the p.d.f explicitly dependent on the choice of
302/// observables. This method instructs RooAddPdf to freeze the
303/// interpretation of the coefficients to be done in the given set of
304/// observables. If frozen, fractions are automatically transformed
305/// from the reference normalization set to the contextual normalization
306/// set by ratios of integrals.
307
309{
310 if (refCoefNorm.empty()) {
311 return ;
312 }
313
315 _refCoefNorm.add(refCoefNorm) ;
316
318}
319
320
321
322////////////////////////////////////////////////////////////////////////////////
323/// By default, fraction coefficients are assumed to refer to the default
324/// fit range. This makes the shape of a RooAddPdf
325/// explicitly dependent on the range of the observables. Calling this function
326/// allows for a range-independent definition of the fractions, because it
327/// ties all coefficients to the given
328/// named range. If the normalisation range is different
329/// from this reference range, the appropriate fraction coefficients
330/// are automatically calculated from the reference fractions by
331/// integrating over the ranges, and comparing these integrals.
332
333void RooAddPdf::fixCoefRange(const char* rangeName)
334{
335 auto* newNamePtr = const_cast<TNamed*>(RooNameReg::ptr(rangeName));
336 if(newNamePtr != _refCoefRangeName) {
338 }
339 _refCoefRangeName = newNamePtr;
340}
341
342
343
344////////////////////////////////////////////////////////////////////////////////
345/// Retrieve cache element for the computation of the PDF normalisation.
346/// \param[in] nset Current normalisation set (integration over these variables yields 1).
347/// \param[in] iset Integration set. Variables to be integrated over (if integrations are performed).
348/// \param[in] rangeName Reference range for the integrals.
349///
350/// If a cache element does not exist, create and fill it on the fly. The cache also contains
351/// - Supplemental normalization terms (in case not all added p.d.f.s have the same observables)
352/// - Projection integrals to calculate transformed fraction coefficients when a frozen reference frame is provided
353/// - Projection integrals for similar transformations when a frozen reference range is provided.
354
356{
357 // Check if cache already exists
358 auto cache = static_cast<AddCacheElem*>(_projCacheMgr.getObj(nset,iset,0,normRange()));
359 if (cache) {
360 return cache ;
361 }
362
363 //Create new cache
364 cache = new AddCacheElem{*this, _pdfList, _coefList, nset, iset, _refCoefNorm,
367 //std::cout << std::endl;
368 //cache->print();
369 //std::cout << std::endl;
370
371 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(normRange())) ;
372
373 return cache;
374}
375
376
377////////////////////////////////////////////////////////////////////////////////
378/// Update the coefficient values in the given cache element: calculate new remainder
379/// fraction, normalize fractions obtained from extended ML terms to unity, and
380/// multiply the various range and dimensional corrections needed in the
381/// current use context.
382///
383/// param[in] cache The cache element for the given normalization set that
384/// stores the supplementary normalization values and
385/// projection-related objects.
386/// param[in] nset The set of variables to normalize over.
387/// param[in] syncCoefValues If the initial values of the coefficients still
388/// need to be copied from the `_coefList` elements to
389/// the `_coefCache`. True by default.
390
391void RooAddPdf::updateCoefficients(AddCacheElem& cache, const RooArgSet* nset, bool syncCoefValues) const
392{
393 _coefCache.resize(_pdfList.size());
394 if(syncCoefValues) {
395 for(std::size_t i = 0; i < _coefList.size(); ++i) {
396 _coefCache[i] = static_cast<RooAbsReal const&>(_coefList[i]).getVal(nset);
397 }
398 }
401}
402
403////////////////////////////////////////////////////////////////////////////////
404/// Look up projection cache and per-PDF norm sets. If a PDF doesn't have a special
405/// norm set, use the `defaultNorm`. If `defaultNorm == nullptr`, use the member
406/// _normSet.
407std::pair<const RooArgSet*, AddCacheElem*> RooAddPdf::getNormAndCache(const RooArgSet* nset) const {
408
409 // Treat empty normalization set and nullptr the same way.
410 if(nset && nset->empty()) nset = nullptr;
411
412 if (nset == nullptr) {
413 if (!_refCoefNorm.empty()) {
414 nset = &_refCoefNorm ;
415 }
416 }
417
418 // A RooAddPdf needs to have a normalization set defined, otherwise its
419 // coefficient will not be uniquely defined. Its shape depends on the
420 // normalization provided. Un-normalized calls to RooAddPdf can happen in
421 // Roofit, when printing the pdf's or when computing integrals. In these case,
422 // if the pdf has a normalization set previously defined (i.e. stored as a
423 // datamember in _copyOfLastNormSet) it should use it by default when the pdf
424 // is evaluated without passing a normalizations set (in pdf->getVal(nullptr) )
425 // In the case of no pre-defined normalization set exists, a warning will be
426 // produced, since the obtained value will be arbitrary. Note that to avoid
427 // unnecessary warning messages, when calling RooAbsPdf::printValue or
428 // RooAbsPdf::graphVizTree, the printing of the warning messages for the
429 // RooFit::Eval topic is explicitly disabled.
430 {
431 // If nset is still nullptr, get the pointer to a copy of the last-used
432 // normalization set. It nset is not nullptr, check whether the copy of
433 // the last-used normalization set needs an update.
434 if(nset == nullptr) {
435 nset = _copyOfLastNormSet.get();
436 } else if(nset->uniqueId() != _idOfLastUsedNormSet) {
437 _copyOfLastNormSet = std::make_unique<const RooArgSet>(*nset);
439 }
440
441 // If nset is STILL nullptr, print a warning.
442 if (nset == nullptr) {
443 coutW(Eval) << "Evaluating RooAddPdf " << GetName() << " without a defined normalization set. This can lead to ambiguous "
444 "coefficients definition and incorrect results."
445 << " Use RooAddPdf::fixCoefNormalization(nset) to provide a normalization set for "
446 "defining uniquely RooAddPdf coefficients!"
447 << std::endl;
448 }
449 }
450
451
452 AddCacheElem* cache = getProjCache(nset) ;
453
454 return {nset, cache};
455}
456
457
458////////////////////////////////////////////////////////////////////////////////
459/// Calculate and return the current value
460
461double RooAddPdf::getValV(const RooArgSet* normSet) const
462{
463 auto normAndCache = getNormAndCache(normSet);
464 const RooArgSet* nset = normAndCache.first;
465 AddCacheElem* cache = normAndCache.second;
466 updateCoefficients(*cache, nset);
467
468 // Process change in last data set used
469 bool nsetChanged(false) ;
470 if (!isActiveNormSet(nset) || _norm==0) {
471 nsetChanged = syncNormalization(nset) ;
472 }
473
474 // Do running sum of coef/pdf pairs, calculate lastCoef.
475 if (isValueDirty() || nsetChanged) {
476 _value = 0.0;
477
478 for (unsigned int i=0; i < _pdfList.size(); ++i) {
479 auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
480 double snormVal = 1.;
481 snormVal = cache->suppNormVal(i);
482
483 double pdfVal = pdf.getVal(nset);
484 if (pdf.isSelectedComp()) {
485 _value += pdfVal*_coefCache[i]/snormVal;
486 }
487 }
489 }
490
491 return _value;
492}
493
494
495////////////////////////////////////////////////////////////////////////////////
496/// Compute addition of PDFs in batches.
497void RooAddPdf::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
498{
499 _coefCache.resize(_pdfList.size());
500 for(std::size_t i = 0; i < _coefList.size(); ++i) {
501 auto coefVals = dataMap.at(&_coefList[i]);
502 // We don't support per-event coefficients in this function. If the CPU
503 // mode is used, we can just fall back to the RooAbsReal implementation.
504 // With CUDA, we can't do that because the inputs might be on the device.
505 // That's why we throw an exception then.
506 if(coefVals.size() > 1) {
507 if(stream) {
508 throw std::runtime_error("The RooAddPdf doesn't support per-event coefficients in CUDA mode yet!");
509 }
510 RooAbsReal::computeBatch(stream, output, nEvents, dataMap);
511 return;
512 }
513 _coefCache[i] = coefVals[0];
514 }
515
518 auto normAndCache = getNormAndCache(nullptr);
519 const RooArgSet* nset = normAndCache.first;
520 AddCacheElem* cache = normAndCache.second;
521 // We don't sync the coefficient values from the _coefList to the _coefCache
522 // because we have already done it using the dataMap.
523 updateCoefficients(*cache, nset, /*syncCoefValues=*/false);
524
525 for (unsigned int pdfNo = 0; pdfNo < _pdfList.size(); ++pdfNo)
526 {
527 auto pdf = static_cast<RooAbsPdf*>(&_pdfList[pdfNo]);
528 if (pdf->isSelectedComp())
529 {
530 pdfs.push_back(dataMap.at(pdf));
531 coefs.push_back(_coefCache[pdfNo] / cache->suppNormVal(pdfNo) );
532 }
533 }
535 dispatch->compute(stream, RooBatchCompute::AddPdf, output, nEvents, pdfs, coefs);
536}
537
538
539////////////////////////////////////////////////////////////////////////////////
540/// Reset error counter to given value, limiting the number
541/// of future error messages for this pdf to 'resetValue'
542
544{
546 _coefErrCount = resetValue ;
547}
548
549
550
551////////////////////////////////////////////////////////////////////////////////
552/// Check if PDF is valid for given normalization set.
553/// Coeffient and PDF must be non-overlapping, but pdf-coefficient
554/// pairs may overlap each other
555
557{
559}
560
561
562////////////////////////////////////////////////////////////////////////////////
563/// Determine which part (if any) of given integral can be performed analytically.
564/// If any analytical integration is possible, return integration scenario code
565///
566/// RooAddPdf queries each component PDF for its analytical integration capability of the requested
567/// set ('allVars'). It finds the largest common set of variables that can be integrated
568/// by all components. If such a set exists, it reconfirms that each component is capable of
569/// analytically integrating the common set, and combines the components individual integration
570/// codes into a single integration code valid for RooAddPdf.
571
573 const RooArgSet* normSet, const char* rangeName) const
574{
575
576 RooArgSet allAnalVars(*std::unique_ptr<RooArgSet>{getObservables(allVars)}) ;
577
578 Int_t n(0) ;
579
580 // First iteration, determine what each component can integrate analytically
581 for (const auto pdfArg : _pdfList) {
582 auto pdf = static_cast<const RooAbsPdf *>(pdfArg);
583 RooArgSet subAnalVars ;
584 pdf->getAnalyticalIntegralWN(allVars,subAnalVars,normSet,rangeName) ;
585
586 // Observables that cannot be integrated analytically by this component are dropped from the common list
587 for (const auto arg : allVars) {
588 if (!subAnalVars.find(arg->GetName()) && pdf->dependsOn(*arg)) {
589 allAnalVars.remove(*arg,true,true) ;
590 }
591 }
592 n++ ;
593 }
594
595 // If no observables can be integrated analytically, return code 0 here
596 if (allAnalVars.empty()) {
597 return 0 ;
598 }
599
600
601 // Now retrieve codes for integration over common set of analytically integrable observables for each component
602 n=0 ;
603 std::vector<Int_t> subCode(_pdfList.size());
604 bool allOK(true) ;
605 for (const auto arg : _pdfList) {
606 auto pdf = static_cast<const RooAbsPdf *>(arg);
607 RooArgSet subAnalVars ;
608 auto allAnalVars2 = std::unique_ptr<RooArgSet>{pdf->getObservables(allAnalVars)} ;
609 subCode[n] = pdf->getAnalyticalIntegralWN(*allAnalVars2,subAnalVars,normSet,rangeName) ;
610 if (subCode[n]==0 && !allAnalVars2->empty()) {
611 coutE(InputArguments) << "RooAddPdf::getAnalyticalIntegral(" << GetName() << ") WARNING: component PDF " << pdf->GetName()
612 << " advertises inconsistent set of integrals (e.g. (X,Y) but not X or Y individually."
613 << " Distributed analytical integration disabled. Please fix PDF" << std::endl ;
614 allOK = false ;
615 }
616 n++ ;
617 }
618 if (!allOK) {
619 return 0 ;
620 }
621
622 // Mare all analytically integrated observables as such
623 analVars.add(allAnalVars) ;
624
625 // Store set of variables analytically integrated
626 RooArgSet* intSet = new RooArgSet(allAnalVars) ;
627 Int_t masterCode = _codeReg.store(subCode,intSet)+1 ;
628
629 return masterCode ;
630}
631
632
633
634////////////////////////////////////////////////////////////////////////////////
635/// Return analytical integral defined by given scenario code
636
637double RooAddPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
638{
639 // WVE needs adaptation to handle new rangeName feature
640 if (code==0) {
641 return getVal(normSet) ;
642 }
643
644 // Retrieve analytical integration subCodes and set of observabels integrated over
645 RooArgSet* intSet ;
646 const std::vector<Int_t>& subCode = _codeReg.retrieve(code-1,intSet) ;
647 if (subCode.empty()) {
648 std::stringstream errorMsg;
649 errorMsg << "RooAddPdf::analyticalIntegral(" << GetName() << "): ERROR unrecognized integration code, " << code;
650 coutE(InputArguments) << errorMsg.str() << std::endl;
651 throw std::invalid_argument(errorMsg.str().c_str());
652 }
653
654 cxcoutD(Caching) << "RooAddPdf::aiWN(" << GetName() << ") calling getProjCache with nset = " << (normSet?*normSet:RooArgSet()) << std::endl ;
655
656 if ((normSet==0 || normSet->empty()) && !_refCoefNorm.empty()) {
657// cout << "WVE integration of RooAddPdf without normalization, but have reference set, using ref set for normalization" << std::endl ;
658 normSet = &_refCoefNorm ;
659 }
660
661 AddCacheElem* cache = getProjCache(normSet,intSet);
662 updateCoefficients(*cache,normSet);
663
664 // Calculate the current value of this object
665 double value(0) ;
666
667 // Do running sum of coef/pdf pairs, calculate lastCoef.
668 double snormVal ;
669
670 //cout << "ROP::aIWN updateCoefCache with rangeName = " << (rangeName?rangeName:"<null>") << std::endl ;
671 for (std::size_t i = 0; i < _pdfList.size(); ++i ) {
672 auto pdf = static_cast<const RooAbsPdf*>(_pdfList.at(i));
673
674 if (_coefCache[i]) {
675 snormVal = cache->suppNormVal(i);
676
677 // WVE swap this?
678 double val = pdf->analyticalIntegralWN(subCode[i],normSet,rangeName) ;
679 if (pdf->isSelectedComp()) {
680 value += val*_coefCache[i]/snormVal ;
681 }
682 }
683 }
684
685 return value ;
686}
687
688
689
690////////////////////////////////////////////////////////////////////////////////
691/// Return the number of expected events, which is either the sum of all coefficients
692/// or the sum of the components extended terms, multiplied with the fraction that
693/// is in the current range w.r.t the reference range
694
695double RooAddPdf::expectedEvents(const RooArgSet* nset) const
696{
697 double expectedTotal{0.0};
698
699 cxcoutD(Caching) << "RooAddPdf::expectedEvents(" << GetName() << ") calling getProjCache with nset = " << (nset?*nset:RooArgSet()) << std::endl ;
700 AddCacheElem& cache = *getProjCache(nset) ;
701 updateCoefficients(cache, nset);
702
703 if (cache.doProjection()) {
704
705 for (std::size_t i = 0; i < _pdfList.size(); ++i) {
706 double ncomp = _allExtendable ? static_cast<RooAbsPdf&>(_pdfList[i]).expectedEvents(nset)
707 : static_cast<RooAbsReal&>(_coefList[i]).getVal(nset);
708 expectedTotal += cache.rangeProjScaleFactor(i) * ncomp ;
709
710 }
711
712 } else {
713
714 if (_allExtendable) {
715 for(auto const& arg : _pdfList) {
716 expectedTotal += static_cast<RooAbsPdf*>(arg)->expectedEvents(nset) ;
717 }
718 } else {
719 for(auto const& arg : _coefList) {
720 expectedTotal += static_cast<RooAbsReal*>(arg)->getVal(nset) ;
721 }
722 }
723
724 }
725 return expectedTotal ;
726}
727
728
729
730////////////////////////////////////////////////////////////////////////////////
731/// Interface function used by test statistics to freeze choice of observables
732/// for interpretation of fraction coefficients
733
734void RooAddPdf::selectNormalization(const RooArgSet* depSet, bool force)
735{
736
737 if (!force && !_refCoefNorm.empty()) {
738 return ;
739 }
740
741 if (!depSet) {
743 return ;
744 }
745
746 fixCoefNormalization(*std::unique_ptr<RooArgSet>{getObservables(depSet)}) ;
747}
748
749
750
751////////////////////////////////////////////////////////////////////////////////
752/// Interface function used by test statistics to freeze choice of range
753/// for interpretation of fraction coefficients
754
755void RooAddPdf::selectNormalizationRange(const char* rangeName, bool force)
756{
757 if (!force && _refCoefRangeName) {
758 return ;
759 }
760
761 fixCoefRange(rangeName) ;
762}
763
764
765
766////////////////////////////////////////////////////////////////////////////////
767/// Return specialized context to efficiently generate toy events from RooAddPdfs
768/// return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ; // WVE DEBUG
769
771 const RooArgSet* auxProto, bool verbose) const
772{
773 return RooAddGenContext::create(*this,vars,prototype,auxProto,verbose).release();
774}
775
776
777
778////////////////////////////////////////////////////////////////////////////////
779/// Loop over components for plot sampling hints and merge them if there are multiple
780
781std::list<double>* RooAddPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
782{
783 return RooRealSumPdf::plotSamplingHint(_pdfList, obs, xlo, xhi);
784}
785
786
787////////////////////////////////////////////////////////////////////////////////
788/// Loop over components for plot sampling hints and merge them if there are multiple
789
790std::list<double>* RooAddPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
791{
792 return RooRealSumPdf::binBoundaries(_pdfList, obs, xlo, xhi);
793}
794
795
796////////////////////////////////////////////////////////////////////////////////
797/// If all components that depend on obs are binned, so is their sum.
799{
801}
802
803
804////////////////////////////////////////////////////////////////////////////////
805/// Label OK'ed components of a RooAddPdf with cache-and-track
806
808{
810}
811
812
813
814////////////////////////////////////////////////////////////////////////////////
815/// Customized printing of arguments of a RooAddPdf to more intuitively reflect the contents of the
816/// product operator construction
817
818void RooAddPdf::printMetaArgs(std::ostream& os) const
819{
821}
822
823
824bool RooAddPdf::redirectServersHook(const RooAbsCollection & newServerList, bool mustReplaceAll,
825 bool nameChange, bool isRecursiveStep)
826{
827 // If a server is redirected, the cached normalization set might not point
828 // to the right observables anymore. We need to reset it.
829 _copyOfLastNormSet.reset();
830 return RooAbsPdf::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursiveStep);
831}
#define cxcoutD(a)
Definition: RooMsgService.h:85
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#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
double rangeProjScaleFactor(std::size_t idx) const
Definition: RooAddHelpers.h:40
bool doProjection() const
Definition: RooAddHelpers.h:31
double suppNormVal(std::size_t idx) const
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=nullptr, RooArgSet *set2=nullptr, RooArgSet *set3=nullptr, RooArgSet *set4=nullptr)
Store given arrays of integer codes, and up to four RooArgSets in the registry (each setX pointer may...
void clearValueAndShapeDirty() const
Definition: RooAbsArg.h:596
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:293
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2185
bool isValueDirty() const
Definition: RooAbsArg.h:421
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
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:576
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:673
bool isActiveNormSet(RooArgSet const *normSet) const
Checks if normSet is the currently active normalization set of this PDF, meaning is exactly the same ...
Definition: RooAbsPdf.h:344
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
Int_t _errorCount
Number of errors remaining to print.
Definition: RooAbsPdf.h:380
const char * normRange() const
Definition: RooAbsPdf.h:300
bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursiveStep) override
The cache manager.
Definition: RooAbsPdf.cxx:3603
static Int_t _verboseEval
Definition: RooAbsPdf.h:359
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
virtual void computeBatch(cudaStream_t *, double *output, size_t size, RooFit::Detail::DataMap const &) const
Base function for computing multiple values of a RooAbsReal.
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...
static void updateCoefficients(RooAbsPdf const &addPdf, std::vector< double > &coefCache, RooArgList const &pdfList, bool haveLastCoef, AddCacheElem &cache, const RooArgSet *nset, RooArgSet const &refCoefNormSet, bool allExtendable, int &coefErrCount)
Update the RooAddPdf coefficients for a given normalization set and projection configuration.
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:128
bool _allExtendable
Flag indicating if all PDF components are extendable.
Definition: RooAddPdf.h:132
RooAICRegistry _codeReg
! Registry of component analytical integration codes
Definition: RooAddPdf.h:125
RooFit::UniqueId< RooArgSet >::Value_t _idOfLastUsedNormSet
!
Definition: RooAddPdf.h:141
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:637
std::unique_ptr< const RooArgSet > _copyOfLastNormSet
!
Definition: RooAddPdf.h:142
void updateCoefficients(AddCacheElem &cache, const RooArgSet *nset, bool syncCoefValues=true) const
Update the coefficient values in the given cache element: calculate new remainder fraction,...
Definition: RooAddPdf.cxx:391
Int_t _coefErrCount
! Coefficient error counter
Definition: RooAddPdf.h:135
bool _haveLastCoef
Flag indicating if last PDFs coefficient was supplied in the ctor.
Definition: RooAddPdf.h:131
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:734
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:818
void finalizeConstruction()
Definition: RooAddPdf.cxx:102
void setCacheAndTrackHints(RooArgSet &) override
Label OK'ed components of a RooAddPdf with cache-and-track.
Definition: RooAddPdf.cxx:807
bool _recursive
Flag indicating is fractions are treated recursively.
Definition: RooAddPdf.h:133
RooObjCacheManager _projCacheMgr
Definition: RooAddPdf.h:106
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:770
bool checkObservables(const RooArgSet *nset) const override
Check if PDF is valid for given normalization set.
Definition: RooAddPdf.cxx:556
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:308
std::pair< const RooArgSet *, AddCacheElem * > getNormAndCache(const RooArgSet *nset) const
Look up projection cache and per-PDF norm sets.
Definition: RooAddPdf.cxx:407
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:755
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:572
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:543
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:695
double getValV(const RooArgSet *set=nullptr) const override
Calculate and return the current value.
Definition: RooAddPdf.cxx:461
void fixCoefRange(const char *rangeName)
By default, fraction coefficients are assumed to refer to the default fit range.
Definition: RooAddPdf.cxx:333
AddCacheElem * getProjCache(const RooArgSet *nset, const RooArgSet *iset=nullptr) const
Manager of cache with coefficient projections and transformations.
Definition: RooAddPdf.cxx:355
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:781
RooListProxy _pdfList
List of component PDFs.
Definition: RooAddPdf.h:127
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:497
bool isBinnedDistribution(const RooArgSet &obs) const override
If all components that depend on obs are binned, so is their sum.
Definition: RooAddPdf.cxx:798
bool redirectServersHook(const RooAbsCollection &, bool, bool, bool) override
The cache manager.
Definition: RooAddPdf.cxx:824
std::vector< double > _coefCache
! Transient cache with transformed values of coefficients
Definition: RooAddPdf.h:103
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:790
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 const char * str(const TNamed *ptr)
Return C++ string corresponding to given TNamed pointer.
Definition: RooNameReg.h:37
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
RooConstVar & RooConst(double val)
const Int_t n
Definition: legend1.C:16
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:62
static void output()