ROOT   Reference Guide
Searching...
No Matches
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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17//////////////////////////////////////////////////////////////////////////////
19 \ingroup Roofitcore
20
21Efficient 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
67
70#include <RooBatchCompute.h>
71#include <RooDataSet.h>
72#include <RooGenericPdf.h>
73#include <RooGlobalFunc.h>
74#include <RooProduct.h>
75#include <RooRatio.h>
76#include <RooRealConstant.h>
77#include <RooRealProxy.h>
78#include <RooRealSumFunc.h>
79#include <RooRealSumPdf.h>
80#include <RooRealVar.h>
82
84#include "RooFitImplHelpers.h"
85
86#include <ROOT/StringUtils.hxx>
87
88#include <algorithm>
89#include <memory>
90#include <set>
91#include <sstream>
92
94
95
96////////////////////////////////////////////////////////////////////////////////
97/// Dummy constructor
98
100 RooAbsPdf(name,title),
101 _refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,false,false),
102 _projCacheMgr(this,10),
103 _pdfList("!pdfs","List of PDFs",this),
104 _coefList("!coefficients","List of coefficients",this),
105 _coefErrCount{_errorCount}
106{
108}
109
110
112
113 // Two pdfs with the same name are only allowed in the input list if they are
114 // actually the same object.
115 using PdfInfo = std::pair<std::string,RooAbsArg*>;
116 std::set<PdfInfo> seen;
117 for(auto const& pdf : _pdfList) {
118 PdfInfo elem{pdf->GetName(), pdf};
119 auto comp = [&](PdfInfo const& p){ return p.first == elem.first && p.second != elem.second; };
120 auto found = std::find_if(seen.begin(), seen.end(), comp);
121 if(found != seen.end()) {
122 std::stringstream errorMsg;
124 << ") pdf list contains pdfs with duplicate name \"" << pdf->GetName() << "\".";
125 coutE(InputArguments) << errorMsg.str() << std::endl;
126 throw std::invalid_argument(errorMsg.str().c_str());
127 }
128 seen.insert(elem);
129 }
130}
131
132
133////////////////////////////////////////////////////////////////////////////////
134/// Constructor with two PDFs and one coefficient
135
137 RooAbsPdf& pdf1, RooAbsPdf& pdf2, RooAbsReal& coef1) :
139{
143
145}
146
147
148////////////////////////////////////////////////////////////////////////////////
149/// Generic constructor from list of PDFs and list of coefficients.
150/// Each pdf list element (i) is paired with coefficient list element (i).
151/// The number of coefficients must be either equal to the number of PDFs,
152/// in which case extended MLL fitting is enabled, or be one less.
153///
154/// All PDFs must inherit from RooAbsPdf. All coefficients must inherit from RooAbsReal
155///
156/// If the recursiveFraction flag is true, the coefficients are interpreted as recursive
157/// coefficients as explained in the class description.
158
159RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList &inPdfList, const RooArgList &inCoefList,
160 bool recursiveFractions)
162{
163 setRecursiveFraction(recursiveFractions);
164
165 if (inPdfList.size()>inCoefList.size()+1 || inPdfList.size()<inCoefList.size()) {
166 std::stringstream errorMsg;
168 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1.";
169 coutE(InputArguments) << errorMsg.str() << std::endl;
170 throw std::invalid_argument(errorMsg.str().c_str());
171 }
172
173 if (recursiveFractions && inPdfList.size()!=inCoefList.size()+1) {
174 std::stringstream errorMsg;
176 << "): Recursive fractions option can only be used if Npdf=Ncoef+1.";
177 coutE(InputArguments) << errorMsg.str() << std::endl;
178 throw std::invalid_argument(errorMsg.str());
179 }
180
181 // Constructor with N PDFs and N or N-1 coefs
182 RooArgList partinCoefList ;
183
184 auto addRecursiveCoef = [this,&partinCoefList](RooAbsPdf& pdf, RooAbsReal& coef) -> RooAbsReal & {
186 if(partinCoefList.size() == 1) {
187 // The first fraction is the first plain fraction
188 return coef;
189 }
190 // The i-th recursive fraction = (1-f1)*(1-f2)*...(fi) and is calculated from the list (f1,...,fi) by RooRecursiveFraction)
191 std::stringstream rfracName;
192 rfracName << GetName() << "_recursive_fraction_" << pdf.GetName() << "_" << partinCoefList.size();
193 auto rfrac = std::make_unique<RooRecursiveFraction>(rfracName.str().c_str(),"Recursive Fraction",partinCoefList) ;
194 auto & rfracRef = *rfrac;
196 return rfracRef;
197 };
198
199 for (auto i = 0u; i < inCoefList.size(); ++i) {
200 auto coef = dynamic_cast<RooAbsReal*>(inCoefList.at(i));
201 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(i));
202 if (inPdfList.at(i) == nullptr) {
203 std::stringstream errorMsg;
205 << ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1";
206 coutE(InputArguments) << errorMsg.str() << std::endl;
207 throw std::invalid_argument(errorMsg.str());
208 }
209 if (!coef) {
210 std::stringstream errorMsg;
211 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") coefficient " << (coef ? coef->GetName() : "") << " is not of type RooAbsReal, ignored";
212 coutE(InputArguments) << errorMsg.str() << std::endl;
213 throw std::invalid_argument(errorMsg.str());
214 }
215 if (!pdf) {
216 std::stringstream errorMsg;
217 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "") << " is not of type RooAbsPdf, ignored";
218 coutE(InputArguments) << errorMsg.str() << std::endl;
219 throw std::invalid_argument(errorMsg.str());
220 }
222
223 // Process recursive fraction mode separately
225 }
226
227 if (inPdfList.size() == inCoefList.size() + 1) {
228 auto pdf = dynamic_cast<RooAbsPdf*>(inPdfList.at(inCoefList.size()));
229
230 if (!pdf) {
231 coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") last argument " << inPdfList.at(inCoefList.size())->GetName() << " is not of type RooAbsPdf." << std::endl;
232 throw std::invalid_argument("Last argument for RooAddPdf is not a PDF.");
233 }
235
236 // Process recursive fractions mode. Above, we verified that we don't have a last coefficient
237 if (recursiveFractions) {
239 // In recursive mode we always have Ncoef=Npdf, since we added it just above
240 _haveLastCoef=true ;
241 }
242
243 } else {
244 _haveLastCoef=true ;
245 }
246
248}
249
250
251////////////////////////////////////////////////////////////////////////////////
252/// Generic constructor from list of extended PDFs. There are no coefficients as the expected
253/// number of events from each components determine the relative weight of the PDFs.
254///
255/// All PDFs must inherit from RooAbsPdf.
256
259{
260 setAllExtendable(true);
261
262 // Constructor with N PDFs
263 for (const auto pdfArg : inPdfList) {
264 auto pdf = dynamic_cast<const RooAbsPdf*>(pdfArg);
265
266 if (!pdf) {
267 std::stringstream errorMsg;
268 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << (pdf ? pdf->GetName() : "")
269 << " is not of type RooAbsPdf, RooAddPdf constructor call is invalid!";
270 coutE(InputArguments) << errorMsg.str() << std::endl;
271 throw std::invalid_argument(errorMsg.str().c_str());
272 }
273 if (!pdf->canBeExtended()) {
274 std::stringstream errorMsg;
275 errorMsg << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName()
276 << " is not extendable, RooAddPdf constructor call is invalid!";
277 coutE(InputArguments) << errorMsg.str() << std::endl;
278 throw std::invalid_argument(errorMsg.str().c_str());
279 }
281 }
282
284}
285
286
287////////////////////////////////////////////////////////////////////////////////
288/// Copy constructor
289
291 : RooAbsPdf(other, name),
292 _refCoefNorm("!refCoefNorm", this, other._refCoefNorm),
293 _refCoefRangeName((TNamed *)other._refCoefRangeName),
294 _projCacheMgr(other._projCacheMgr, this),
295 _codeReg(other._codeReg),
296 _pdfList("!pdfs", this, other._pdfList),
297 _coefList("!coefficients", this, other._coefList),
298 _haveLastCoef(other._haveLastCoef),
299 _allExtendable(other._allExtendable),
300 _recursive(other._recursive),
301 _coefErrCount(_errorCount)
302{
303
306}
307
308
309////////////////////////////////////////////////////////////////////////////////
310/// By default the interpretation of the fraction coefficients is
311/// performed in the contextual choice of observables. This makes the
312/// shape of the p.d.f explicitly dependent on the choice of
313/// observables. This method instructs RooAddPdf to freeze the
314/// interpretation of the coefficients to be done in the given set of
315/// observables. If frozen, fractions are automatically transformed
316/// from the reference normalization set to the contextual normalization
317/// set by ratios of integrals.
318
320{
321 if (refCoefNorm.empty()) {
322 return ;
323 }
324
325 // Also set an attribute with this information, which is the easiest way to
326 // preserve this in the JSON IO.
327 setStringAttribute("ref_coef_norm", RooHelpers::getColonSeparatedNameString(refCoefNorm, ',').c_str());
328
331
333}
334
336{
338 return _refCoefNorm;
339}
340
341// For the JSON IO, we are not storing the _refCoefNorm directly. Instead, it
342// is stored by names in a string attribute. This function should be called
343// internally before _refCoefNorm is used to materialize it from the attribute
344// if necessary.
346{
347 // _refCoefNorm was already materialized
348 if (!_refCoefNorm.empty())
349 return;
350
351 std::vector<std::string> names;
352 if (auto attrib = getStringAttribute("ref_coef_norm")) {
353 names = ROOT::Split(attrib, ",", /*skipEmpty=*/true);
354 } else {
355 return;
356 }
357
358 RooArgSet refCoefNorm;
359
360 RooArgSet serverSet;
362 for (std::string const &name : names) {
363 if (RooAbsArg *arg = serverSet.find(name.c_str())) {
365 } else {
366 throw std::runtime_error("Internal logic error in RooAddPdf::materializeRefCoefNormFromAttribute()");
367 }
368 }
369
371}
372
373
374////////////////////////////////////////////////////////////////////////////////
375/// By default, fraction coefficients are assumed to refer to the default
376/// fit range. This makes the shape of a RooAddPdf
377/// explicitly dependent on the range of the observables. Calling this function
378/// allows for a range-independent definition of the fractions, because it
379/// ties all coefficients to the given
380/// named range. If the normalisation range is different
381/// from this reference range, the appropriate fraction coefficients
382/// are automatically calculated from the reference fractions by
383/// integrating over the ranges, and comparing these integrals.
384
386{
387 auto* newNamePtr = const_cast<TNamed*>(RooNameReg::ptr(rangeName));
388 if(newNamePtr != _refCoefRangeName) {
390 }
391 _refCoefRangeName = newNamePtr;
392}
393
394
395
396////////////////////////////////////////////////////////////////////////////////
397/// Retrieve cache element for the computation of the PDF normalisation.
398/// \param[in] nset Current normalisation set (integration over these variables yields 1).
399/// \param[in] iset Integration set. Variables to be integrated over (if integrations are performed).
400///
401/// If a cache element does not exist, create and fill it on the fly. The cache also contains
402/// - Supplemental normalization terms (in case not all added p.d.f.s have the same observables)
403/// - Projection integrals to calculate transformed fraction coefficients when a frozen reference frame is provided
404/// - Projection integrals for similar transformations when a frozen reference range is provided.
405
407{
408 // Check if cache already exists
410 if (cache) {
411 return cache ;
412 }
413
414 // Make sure _refCoefNorm is defined
416
417 //Create new cache
418 cache = new AddCacheElem{*this, _pdfList, _coefList, nset, iset, _refCoefNorm,
421
422 _projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(normRange())) ;
423
424 return cache;
425}
426
427
428////////////////////////////////////////////////////////////////////////////////
429/// Update the coefficient values in the given cache element: calculate new remainder
430/// fraction, normalize fractions obtained from extended ML terms to unity, and
431/// multiply the various range and dimensional corrections needed in the
432/// current use context.
433///
434/// param[in] cache The cache element for the given normalization set that
435/// stores the supplementary normalization values and
436/// projection-related objects.
437/// param[in] nset The set of variables to normalize over.
438/// param[in] syncCoefValues If the initial values of the coefficients still
439/// need to be copied from the _coefList elements to
440/// the _coefCache. True by default.
441
443{
444 _coefCache.resize(_pdfList.size());
445 if(syncCoefValues) {
446 for(std::size_t i = 0; i < _coefList.size(); ++i) {
447 _coefCache[i] = static_cast<RooAbsReal const&>(_coefList[i]).getVal(nset);
448 }
449 }
452}
453
454////////////////////////////////////////////////////////////////////////////////
455/// Look up projection cache and per-PDF norm sets. If a PDF doesn't have a special
456/// norm set, use the defaultNorm. If defaultNorm == nullptr, use the member
457/// _normSet.
459
460 // Treat empty normalization set and nullptr the same way.
461 if(nset && nset->empty()) nset = nullptr;
462
463 if (nset == nullptr) {
464 // Make sure _refCoefNorm is defined
466
467 if (!_refCoefNorm.empty()) {
468 nset = &_refCoefNorm ;
469 }
470 }
471
472 // A RooAddPdf needs to have a normalization set defined, otherwise its
473 // coefficient will not be uniquely defined. Its shape depends on the
474 // normalization provided. Un-normalized calls to RooAddPdf can happen in
475 // Roofit, when printing the pdf's or when computing integrals. In these case,
476 // if the pdf has a normalization set previously defined (i.e. stored as a
477 // datamember in _copyOfLastNormSet) it should use it by default when the pdf
478 // is evaluated without passing a normalizations set (in pdf->getVal(nullptr) )
479 // In the case of no pre-defined normalization set exists, a warning will be
480 // produced, since the obtained value will be arbitrary. Note that to avoid
481 // unnecessary warning messages, when calling RooAbsPdf::printValue or
482 // RooAbsPdf::graphVizTree, the printing of the warning messages for the
483 // RooFit::Eval topic is explicitly disabled.
484 {
485 // If nset is still nullptr, get the pointer to a copy of the last-used
486 // normalization set. It nset is not nullptr, check whether the copy of
487 // the last-used normalization set needs an update.
488 if(nset == nullptr) {
489 nset = _copyOfLastNormSet.get();
491 _copyOfLastNormSet = std::make_unique<const RooArgSet>(*nset);
493 }
494
495 // If nset is STILL nullptr, print a warning.
496 if (nset == nullptr) {
497 coutW(Eval) << "Evaluating RooAddPdf " << GetName() << " without a defined normalization set. This can lead to ambiguous "
498 "coefficients definition and incorrect results."
499 << " Use RooAddPdf::fixCoefNormalization(nset) to provide a normalization set for "
501 << std::endl;
502 }
503 }
504
505
506 AddCacheElem* cache = getProjCache(nset) ;
507
508 return {nset, cache};
509}
510
511
512////////////////////////////////////////////////////////////////////////////////
513/// Calculate and return the current value
514
516{
517 auto normAndCache = getNormAndCache(normSet);
518 const RooArgSet* nset = normAndCache.first;
520 updateCoefficients(*cache, nset);
521
522 // Process change in last data set used
523 bool nsetChanged(false) ;
524 if (!isActiveNormSet(nset) || _norm==nullptr) {
525 nsetChanged = syncNormalization(nset) ;
526 }
527
528 // Do running sum of coef/pdf pairs, calculate lastCoef.
529 if (isValueDirty() || nsetChanged) {
530 _value = 0.0;
531
532 for (unsigned int i=0; i < _pdfList.size(); ++i) {
533 auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
534 double snormVal = 1.;
535 snormVal = cache->suppNormVal(i);
536
537 double pdfVal = pdf.getVal(nset);
538 if (pdf.isSelectedComp()) {
539 _value += pdfVal*_coefCache[i]/snormVal;
540 }
541 }
543 }
544
545 return _value;
546}
547
549{
550 ctx.addResult(this, RooRealSumPdf::translateImpl(ctx, this, _pdfList, _coefList, true));
551}
552
553////////////////////////////////////////////////////////////////////////////////
554/// Compute addition of PDFs in batches.
556{
557 std::span<double> output = ctx.output();
558
559 RooBatchCompute::Config config = ctx.config(this);
560
561 _coefCache.resize(_pdfList.size());
562 for(std::size_t i = 0; i < _coefList.size(); ++i) {
563 auto coefVals = ctx.at(&_coefList[i]);
564 // We don't support per-event coefficients in this function. If the CPU
565 // mode is used, we can just fall back to the RooAbsReal implementation.
566 // With CUDA, we can't do that because the inputs might be on the device.
567 // That's why we throw an exception then.
568 if(coefVals.size() > 1) {
569 if (config.useCuda()) {
570 throw std::runtime_error("The RooAddPdf doesn't support per-event coefficients in CUDA mode yet!");
571 }
573 return;
574 }
575 _coefCache[i] = coefVals[0];
576 }
577
578 std::vector<std::span<const double>> pdfs;
579 std::vector<double> coefs;
581 // We don't sync the coefficient values from the _coefList to the _coefCache
582 // because we have already done it using the ctx.
583 updateCoefficients(*cache, nullptr, /*syncCoefValues=*/false);
584
585 for (unsigned int pdfNo = 0; pdfNo < _pdfList.size(); ++pdfNo)
586 {
587 auto pdf = static_cast<RooAbsPdf*>(&_pdfList[pdfNo]);
588 if (pdf->isSelectedComp())
589 {
590 pdfs.push_back(ctx.at(pdf));
591 coefs.push_back(_coefCache[pdfNo] / cache->suppNormVal(pdfNo) );
592 }
593 }
595}
596
597
598////////////////////////////////////////////////////////////////////////////////
599/// Reset error counter to given value, limiting the number
600/// of future error messages for this pdf to 'resetValue'
601
603{
605 _coefErrCount = resetValue ;
606}
607
608
609
610////////////////////////////////////////////////////////////////////////////////
611/// Check if PDF is valid for given normalization set.
612/// Coefficient and PDF must be non-overlapping, but pdf-coefficient
613/// pairs may overlap each other
614
616{
618}
619
620
621////////////////////////////////////////////////////////////////////////////////
622/// Determine which part (if any) of given integral can be performed analytically.
623/// If any analytical integration is possible, return integration scenario code
624///
625/// RooAddPdf queries each component PDF for its analytical integration capability of the requested
626/// set ('allVars'). It finds the largest common set of variables that can be integrated
627/// by all components. If such a set exists, it reconfirms that each component is capable of
628/// analytically integrating the common set, and combines the components individual integration
629/// codes into a single integration code valid for RooAddPdf.
630
632 const RooArgSet* normSet, const char* rangeName) const
633{
634 // Make sure _refCoefNorm is defined
636
637 RooArgSet allAnalVars(*std::unique_ptr<RooArgSet>{getObservables(allVars)}) ;
638
639 Int_t n(0) ;
640
641 // First iteration, determine what each component can integrate analytically
642 for (const auto pdfArg : _pdfList) {
643 auto pdf = static_cast<const RooAbsPdf *>(pdfArg);
644 RooArgSet subAnalVars ;
645 pdf->getAnalyticalIntegralWN(allVars,subAnalVars,normSet,rangeName) ;
646
647 // Observables that cannot be integrated analytically by this component are dropped from the common list
648 for (const auto arg : allVars) {
649 if (!subAnalVars.find(arg->GetName()) && pdf->dependsOn(*arg)) {
650 allAnalVars.remove(*arg,true,true) ;
651 }
652 }
653 n++ ;
654 }
655
656 // If no observables can be integrated analytically, return code 0 here
657 if (allAnalVars.empty()) {
658 return 0 ;
659 }
660
661
662 // Now retrieve codes for integration over common set of analytically integrable observables for each component
663 n=0 ;
664 std::vector<Int_t> subCode(_pdfList.size());
665 bool allOK(true) ;
666 for (const auto arg : _pdfList) {
667 auto pdf = static_cast<const RooAbsPdf *>(arg);
668 RooArgSet subAnalVars ;
669 auto allAnalVars2 = std::unique_ptr<RooArgSet>{pdf->getObservables(allAnalVars)} ;
670 subCode[n] = pdf->getAnalyticalIntegralWN(*allAnalVars2,subAnalVars,normSet,rangeName) ;
671 if (subCode[n]==0 && !allAnalVars2->empty()) {
672 coutE(InputArguments) << "RooAddPdf::getAnalyticalIntegral(" << GetName() << ") WARNING: component PDF " << pdf->GetName()
673 << " advertises inconsistent set of integrals (e.g. (X,Y) but not X or Y individually."
674 << " Distributed analytical integration disabled. Please fix PDF" << std::endl ;
675 allOK = false ;
676 }
677 n++ ;
678 }
679 if (!allOK) {
680 return 0 ;
681 }
682
683 // Mare all analytically integrated observables as such
685
686 // Store set of variables analytically integrated
687 RooArgSet* intSet = new RooArgSet(allAnalVars) ;
688 Int_t masterCode = _codeReg.store(subCode,intSet)+1 ;
689
690 return masterCode ;
691}
692
693
694
695////////////////////////////////////////////////////////////////////////////////
696/// Return analytical integral defined by given scenario code
697
698double RooAddPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
699{
700 // WVE needs adaptation to handle new rangeName feature
701 if (code==0) {
702 return getVal(normSet) ;
703 }
704
705 // Retrieve analytical integration subCodes and set of observabels integrated over
706 RooArgSet* intSet ;
707 const std::vector<Int_t>& subCode = _codeReg.retrieve(code-1,intSet) ;
708 if (subCode.empty()) {
709 std::stringstream errorMsg;
710 errorMsg << "RooAddPdf::analyticalIntegral(" << GetName() << "): ERROR unrecognized integration code, " << code;
711 coutE(InputArguments) << errorMsg.str() << std::endl;
712 throw std::invalid_argument(errorMsg.str().c_str());
713 }
714
715 cxcoutD(Caching) << "RooAddPdf::aiWN(" << GetName() << ") calling getProjCache with nset = " << (normSet?*normSet:RooArgSet()) << std::endl ;
716
717 if ((normSet==nullptr || normSet->empty()) && !_refCoefNorm.empty()) {
718// cout << "WVE integration of RooAddPdf without normalization, but have reference set, using ref set for normalization" << std::endl ;
719 normSet = &_refCoefNorm ;
720 }
721
723 updateCoefficients(*cache,normSet);
724
725 // Calculate the current value of this object
726 double value(0) ;
727
728 // Do running sum of coef/pdf pairs, calculate lastCoef.
729 double snormVal ;
730
731 //cout << "ROP::aIWN updateCoefCache with rangeName = " << (rangeName?rangeName:"<null>") << std::endl ;
732 for (std::size_t i = 0; i < _pdfList.size(); ++i ) {
733 auto pdf = static_cast<const RooAbsPdf*>(_pdfList.at(i));
734
735 if (_coefCache[i]) {
736 snormVal = cache->suppNormVal(i);
737
738 // WVE swap this?
739 double val = pdf->analyticalIntegralWN(subCode[i],normSet,rangeName) ;
740 if (pdf->isSelectedComp()) {
741 value += val*_coefCache[i]/snormVal ;
742 }
743 }
744 }
745
746 return value ;
747}
748
749
750
751////////////////////////////////////////////////////////////////////////////////
752/// Return the number of expected events, which is either the sum of all coefficients
753/// or the sum of the components extended terms, multiplied with the fraction that
754/// is in the current range w.r.t the reference range
755
757{
758 double expectedTotal{0.0};
759
760 cxcoutD(Caching) << "RooAddPdf::expectedEvents(" << GetName() << ") calling getProjCache with nset = " << (nset?*nset:RooArgSet()) << std::endl ;
761 AddCacheElem& cache = *getProjCache(nset) ;
762 updateCoefficients(cache, nset);
763
764 if (cache.doProjection()) {
765
766 for (std::size_t i = 0; i < _pdfList.size(); ++i) {
767 double ncomp = _allExtendable ? static_cast<RooAbsPdf&>(_pdfList[i]).expectedEvents(nset)
768 : static_cast<RooAbsReal&>(_coefList[i]).getVal(nset);
769 expectedTotal += cache.rangeProjScaleFactor(i) * ncomp ;
770
771 }
772
773 } else {
774
775 if (_allExtendable) {
776 for(auto const& arg : _pdfList) {
777 expectedTotal += static_cast<RooAbsPdf*>(arg)->expectedEvents(nset) ;
778 }
779 } else {
780 for(auto const& arg : _coefList) {
781 expectedTotal += static_cast<RooAbsReal*>(arg)->getVal(nset) ;
782 }
783 }
784
785 }
786 return expectedTotal ;
787}
788
789
791{
792 std::unique_ptr<RooAbsReal> out;
793
794 auto name = std::string(GetName()) + "_expectedEvents";
795 if (_allExtendable) {
796 RooArgSet sumSet;
797 for (auto *pdf : static_range_cast<RooAbsPdf *>(_pdfList)) {
799 }
800 out = std::make_unique<RooAddition>(name.c_str(), name.c_str(), sumSet);
802 } else {
803 out = std::make_unique<RooAddition>(name.c_str(), name.c_str(), _coefList);
804 }
805
806 RooArgList prodList;
807
808 // Make sure _refCoefNorm is defined
810
811 if (!_allExtendable) {
812 // If the _refCoefNorm is empty or it's equal to normSet anyway, this is not
813 // a conditional pdf and we don't need to do any transformation. See also
814 // RooAddPdf::compileForNormSet() for more explanations on a similar logic.
815 if (!_refCoefNorm.empty() && !nset->equals(_refCoefNorm)) {
817 }
818
819 // Optionally multiply with fractional normalization. I this case, we
820 // replace the original factor stored in "out".
821 if (!_normRange.IsNull()) {
822 std::unique_ptr<RooAbsReal> owner;
823 RooArgList terms;
824 // The integrals own each other in a chain. We do this because it's
825 // not possible to add two objects with the same name via
826 // addOwnedComponents(), and it happens in some user models that some
827 // component pdfs are the same. Hence, the integrals might share names
828 // too and we can't add them all in one go as owned objects of the
829 // final integral sum.
830 for (auto *pdf : static_range_cast<RooAbsPdf *>(_pdfList)) {
831 auto next = std::unique_ptr<RooAbsReal>{pdf->createIntegral(*nset, *nset, _normRange)};
833 if (owner)
835 owner = std::move(next);
836 }
837 auto fracIntegName = std::string(GetName()) + "_integSum";
838 auto fracInteg =
839 std::make_unique<RooRealSumFunc>(fracIntegName.c_str(), fracIntegName.c_str(), _coefList, terms);
841
842 out = std::move(fracInteg);
843 }
844 }
845
846 std::string finalName = std::string(out->GetName()) + "_finalized";
847 if (prodList.empty()) {
848 // If there are no additional factors, just return the single factor we have
849 return out;
850 } else {
852 }
853 auto finalOut = std::make_unique<RooProduct>(finalName.c_str(), finalName.c_str(), prodList);
855 return finalOut;
856}
857
858
859////////////////////////////////////////////////////////////////////////////////
860/// Interface function used by test statistics to freeze choice of observables
861/// for interpretation of fraction coefficients
862
863void RooAddPdf::selectNormalization(const RooArgSet* depSet, bool force)
864{
865 // Make sure _refCoefNorm is defined
867
868 if (!force && !_refCoefNorm.empty()) {
869 return ;
870 }
871
872 if (!depSet) {
874 return ;
875 }
876
877 fixCoefNormalization(*std::unique_ptr<RooArgSet>{getObservables(depSet)}) ;
878}
879
880
881
882////////////////////////////////////////////////////////////////////////////////
883/// Interface function used by test statistics to freeze choice of range
884/// for interpretation of fraction coefficients
885
886void RooAddPdf::selectNormalizationRange(const char* rangeName, bool force)
887{
888 if (!force && _refCoefRangeName) {
889 return ;
890 }
891
892 fixCoefRange(rangeName) ;
893}
894
895
896
897////////////////////////////////////////////////////////////////////////////////
898/// Return specialized context to efficiently generate toy events from RooAddPdfs
899/// return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ; // WVE DEBUG
900
902 const RooArgSet* auxProto, bool verbose) const
903{
905}
906
907
908
909////////////////////////////////////////////////////////////////////////////////
910/// Loop over components for plot sampling hints and merge them if there are multiple
911
912std::list<double>* RooAddPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
913{
914 return RooRealSumPdf::plotSamplingHint(_pdfList, obs, xlo, xhi);
915}
916
917
918////////////////////////////////////////////////////////////////////////////////
919/// Loop over components for plot sampling hints and merge them if there are multiple
920
921std::list<double>* RooAddPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
922{
923 return RooRealSumPdf::binBoundaries(_pdfList, obs, xlo, xhi);
924}
925
926
927////////////////////////////////////////////////////////////////////////////////
928/// If all components that depend on obs are binned, so is their sum.
930{
932}
933
934
935////////////////////////////////////////////////////////////////////////////////
936/// Label OK'ed components of a RooAddPdf with cache-and-track
937
939{
941}
942
943
944
945////////////////////////////////////////////////////////////////////////////////
946/// Customized printing of arguments of a RooAddPdf to more intuitively reflect the contents of the
947/// product operator construction
948
950{
952}
953
954
955bool RooAddPdf::redirectServersHook(const RooAbsCollection & newServerList, bool mustReplaceAll,
956 bool nameChange, bool isRecursiveStep)
957{
958 // If a server is redirected, the cached normalization set might not point
959 // to the right observables anymore. We need to reset it.
960 _copyOfLastNormSet.reset();
961 return RooAbsPdf::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursiveStep);
962}
963
964
965std::unique_ptr<RooAbsArg>
967{
968 // Make sure _refCoefNorm is defined
970
971 auto newArg = std::unique_ptr<RooAbsReal>{static_cast<RooAbsReal *>(Clone())};
972 ctx.markAsCompiled(*newArg);
973
974 // In case conditional observables, e.g. p(x|y), the _refCoefNorm is set to
975 // all observables (x, y) and the normSet doesn't contain the conditional
976 // observables (so it only contains x in this example).
977
978 // If the _refCoefNorm is empty or it's equal to normSet anyway, this is not
979 // a conditional pdf and we don't need to do any transformation.
980 if(_refCoefNorm.empty() || normSet.equals(_refCoefNorm)) {
981 ctx.compileServers(*newArg, normSet);
982 return newArg;
983 }
984
985 // In the conditional case, things become more complicated. The original
986 // getValV() method is covering this case with very complicated logic,
987 // caching multiple new RooFit objects to scale the individual coefficients
989 //
990 // However, it's not complicated what we need to do mathematically:
991 //
992 // Since:
993 // 1. p(x, y) = p(x | y) * p(y)
994 // 2. p(y) = Integral of p(x, y) over x
995 //
996 // We conclude:
997 // p(x, y)
998 // p(x | y) = --------------------------
999 // Integral of p(x, y) over x
1000 //
1001 // What follows is the implementation of this formula in RooFit. By doing
1002 // this here in compileForNormSet(), we don't invoke the old RooAddPdf
1003 // projection caches (note that no conditional pdfs are on the right hand
1004 // side of the equation).
1005 std::string finalName = std::string(GetName()) + "_conditional";
1006 std::unique_ptr<RooAbsReal> denom{newArg->createIntegral(normSet, _refCoefNorm)};
1007 auto finalArg = std::make_unique<RooGenericPdf>(finalName.c_str(), "@0/@1", RooArgList{*newArg, *denom});
1008 ctx.compileServers(*denom, _refCoefNorm);
1009 ctx.markAsCompiled(*denom);
1010 ctx.markAsCompiled(*finalArg);
1011 ctx.compileServers(*newArg, _refCoefNorm);
1014 return finalArg;
1015}
#define cxcoutD(a)
#define coutW(a)
#define coutE(a)
#define TRACE_CREATE
Definition RooTrace.h:23
#define ClassImp(name)
Definition Rtypes.h:377
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
bool doProjection() const
double suppNormVal(std::size_t idx) const
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...
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
void clearValueAndShapeDirty() const
Definition RooAbsArg.h:562
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.
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
Take ownership of the contents of 'comps'.
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
bool isValueDirty() const
Definition RooAbsArg.h:384
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition RooAbsArg.h:89
Abstract container object that can hold multiple RooAbsArg objects.
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
RooFit::UniqueId< RooAbsCollection > const & uniqueId() const
Returns a unique ID that is different for every instantiated RooAbsCollection.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
Storage_t const & get() const
const char * GetName() const override
Returns name of object.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for generator contexts of RooAbsPdf objects.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
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...
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...
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:300
TString _normRange
Normalization range.
Definition RooAbsPdf.h:343
RooAbsReal * _norm
Definition RooAbsPdf.h:320
const char * normRange() const
Definition RooAbsPdf.h:251
bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursiveStep) override
The cache manager.
static Int_t _verboseEval
Definition RooAbsPdf.h:315
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
double _value
Cache for current value of object.
Definition RooAbsReal.h:537
virtual void doEval(RooFit::EvalContext &) 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.
Efficient implementation of a sum of PDFs of the form.
RooListProxy _coefList
List of coefficients.
bool _allExtendable
Flag indicating if all PDF components are extendable.
void doEval(RooFit::EvalContext &) const override
Compute addition of PDFs in batches.
RooAICRegistry _codeReg
! Registry of component analytical integration codes
RooFit::UniqueId< RooArgSet >::Value_t _idOfLastUsedNormSet
!
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Return analytical integral defined by given scenario code.
std::unique_ptr< const RooArgSet > _copyOfLastNormSet
!
void updateCoefficients(AddCacheElem &cache, const RooArgSet *nset, bool syncCoefValues=true) const
Update the coefficient values in the given cache element: calculate new remainder fraction,...
void translate(RooFit::Detail::CodeSquashContext &ctx) const override
This function defines a translation for each RooAbsReal based object that can be used to express the ...
Int_t _coefErrCount
! Coefficient error counter
void setRecursiveFraction(bool recursiveFraction)
bool _haveLastCoef
Flag indicating if last PDFs coefficient was supplied in the constructor.
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...
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooAddPdf to more intuitively reflect the contents of the produ...
void finalizeConstruction()
void setCacheAndTrackHints(RooArgSet &) override
Label OK'ed components of a RooAddPdf with cache-and-track.
RooObjCacheManager _projCacheMgr
void materializeRefCoefNormFromAttribute() const
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...
bool checkObservables(const RooArgSet *nset) const override
Check if PDF is valid for given normalization set.
void fixCoefNormalization(const RooArgSet &refCoefNorm)
By default the interpretation of the fraction coefficients is performed in the contextual choice of o...
std::pair< const RooArgSet *, AddCacheElem * > getNormAndCache(const RooArgSet *nset) const
Look up projection cache and per-PDF norm sets.
RooSetProxy _refCoefNorm
Reference observable set for coefficient interpretation.
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...
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.
void setAllExtendable(bool allExtendable)
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...
double expectedEvents(const RooArgSet *nset) const override
Return expected number of events for extended likelihood calculation, which is the sum of all coeffic...
double getValV(const RooArgSet *set=nullptr) const override
Calculate and return the current value.
std::unique_ptr< RooAbsArg > compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const override
void fixCoefRange(const char *rangeName)
By default, fraction coefficients are assumed to refer to the default fit range.
AddCacheElem * getProjCache(const RooArgSet *nset, const RooArgSet *iset=nullptr) const
Manager of cache with coefficient projections and transformations.
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.
RooListProxy _pdfList
List of component PDFs.
TNamed * _refCoefRangeName
Reference range name for coefficient interpretation.
std::unique_ptr< RooAbsReal > createExpectedEventsFunc(const RooArgSet *nset) const override
Returns an object that represents the expected number of events for a given normalization set,...
bool isBinnedDistribution(const RooArgSet &obs) const override
If all components that depend on obs are binned, so is their sum.
bool redirectServersHook(const RooAbsCollection &, bool, bool, bool) override
The cache manager.
std::vector< double > _coefCache
! Transient cache with transformed values of coefficients
const RooArgSet & getCoefNormalization() const
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Loop over components for plot sampling hints and merge them if there are multiple.
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:24
Minimal configuration struct to steer the evaluation of a single node with the RooBatchCompute librar...
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...
Container class to hold unbinned data.
Definition RooDataSet.h:33
A class to maintain the context for squashing of RooFit models into code.
void addResult(RooAbsArg const *key, std::string const &value)
A function to save an expression that includes/depends on the result of the input node.
void markAsCompiled(RooAbsArg &arg) const
void compileServers(RooAbsArg &arg, RooArgSet const &normSet)
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
static const char * str(const TNamed *ptr)
Return C++ string corresponding to given TNamed pointer.
Definition RooNameReg.h:39
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
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...
static std::string translateImpl(RooFit::Detail::CodeSquashContext &ctx, RooAbsArg const *klass, RooArgList const &funcList, RooArgList const &coefList, bool normalize=false)
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
Bool_t IsNull() const
Definition TString.h:414
RooConstVar & RooConst(double val)
const Int_t n
Definition legend1.C:16
std::vector< std::string > Split(std::string_view str, std::string_view delims, bool skipEmpty=false)
Splits a string at each character in delims.
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})
void getSortedComputationGraph(RooAbsArg const &func, RooArgSet &out)
std::string getColonSeparatedNameString(RooArgSet const &argSet, char delim=':')
unsigned long Value_t
Definition UniqueId.h:41
static void output()