Logo ROOT  
Reference Guide
RooProdPdf.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\file RooProdPdf.cxx
19\class RooProdPdf
20\ingroup Roofitcore
21
22RooProdPdf is an efficient implementation of a product of PDFs of the form
23\f[ \prod_{i=1}^{N} \mathrm{PDF}_i (x, \ldots) \f]
24
25PDFs may share observables. If that is the case any irreducible subset
26of PDFs that share observables will be normalised with explicit numeric
27integration as any built-in normalisation will no longer be valid.
28
29Alternatively, products using conditional PDFs can be defined, *e.g.*
30
31\f[ F(x|y) \cdot G(y), \f]
32
33meaning a PDF \f$ F(x) \f$ **given** \f$ y \f$ and a PDF \f$ G(y) \f$.
34In this construction, \f$ F \f$ is only
35normalised w.r.t \f$ x\f$, and \f$ G \f$ is normalised w.r.t \f$ y \f$. The product in this construction
36is properly normalised.
37
38If exactly one of the component PDFs supports extended likelihood fits, the
39product will also be usable in extended mode, returning the number of expected
40events from the extendable component PDF. The extendable component does not
41have to appear in any specific place in the list.
42**/
43
44#include "RooProdPdf.h"
45#include "RooBatchCompute.h"
46#include "RooRealProxy.h"
47#include "RooProdGenContext.h"
48#include "RooGenProdProj.h"
49#include "RooProduct.h"
50#include "RooNameReg.h"
51#include "RooMsgService.h"
52#include "RooFormulaVar.h"
53#include "RooRealVar.h"
54#include "RooAddition.h"
55#include "RooGlobalFunc.h"
56#include "RooConstVar.h"
57#include "RooWorkspace.h"
58#include "RooRangeBoolean.h"
59#include "RooCustomizer.h"
60#include "RooRealIntegral.h"
61#include "RooTrace.h"
62#include "RooHelpers.h"
63#include "RooBatchCompute.h"
64#include "strtok.h"
65
66#include <cstring>
67#include <sstream>
68#include <algorithm>
69
70#ifndef _WIN32
71#include <strings.h>
72#endif
73
74using namespace std;
75
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// Default constructor
81
83 _cacheMgr(this,10),
84 _cutOff(0),
85 _extendedIndex(-1),
86 _useDefaultGen(false),
87 _refRangeName(0),
88 _selfNorm(true)
89{
90 // Default constructor
92}
93
94
95
96////////////////////////////////////////////////////////////////////////////////
97/// Dummy constructor
98
99RooProdPdf::RooProdPdf(const char *name, const char *title, double cutOff) :
100 RooAbsPdf(name,title),
101 _cacheMgr(this,10),
102 _genCode(10),
103 _cutOff(cutOff),
104 _pdfList("!pdfs","List of PDFs",this),
105 _extendedIndex(-1),
106 _useDefaultGen(false),
107 _refRangeName(0),
108 _selfNorm(true)
109{
111}
112
113
114
115////////////////////////////////////////////////////////////////////////////////
116/// Constructor with 2 PDFs (most frequent use case).
117///
118/// The optional cutOff parameter can be used as a speed optimization if
119/// one or more of the PDF have sizable regions with very small values,
120/// which would pull the entire product of PDFs to zero in those regions.
121///
122/// After each PDF multiplication, the running product is compared with
123/// the cutOff parameter. If the running product is smaller than the
124/// cutOff value, the product series is terminated and remaining PDFs
125/// are not evaluated.
126///
127/// There is no magic value of the cutOff, the user should experiment
128/// to find the appropriate balance between speed and precision.
129/// If a cutoff is specified, the PDFs most likely to be small should
130/// be put first in the product. The default cutOff value is zero.
131///
132
133RooProdPdf::RooProdPdf(const char *name, const char *title,
134 RooAbsPdf& pdf1, RooAbsPdf& pdf2, double cutOff) :
135 RooAbsPdf(name,title),
136 _cacheMgr(this,10),
137 _genCode(10),
138 _cutOff(cutOff),
139 _pdfList("!pdfs","List of PDFs",this),
140 _extendedIndex(-1),
141 _useDefaultGen(false),
142 _refRangeName(0),
143 _selfNorm(true)
144{
145 _pdfList.add(pdf1) ;
146 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>("nset")) ;
147 if (pdf1.canBeExtended()) {
148 _extendedIndex = _pdfList.index(&pdf1) ;
149 }
150
151 _pdfList.add(pdf2) ;
152 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>("nset")) ;
153
154 if (pdf2.canBeExtended()) {
155 if (_extendedIndex>=0) {
156 // Protect against multiple extended terms
157 coutW(InputArguments) << "RooProdPdf::RooProdPdf(" << GetName()
158 << ") multiple components with extended terms detected,"
159 << " product will not be extendible." << endl ;
160 _extendedIndex=-1 ;
161 } else {
163 }
164 }
166}
167
168
169
170////////////////////////////////////////////////////////////////////////////////
171/// Constructor from a list of PDFs.
172///
173/// The optional cutOff parameter can be used as a speed optimization if
174/// one or more of the PDF have sizable regions with very small values,
175/// which would pull the entire product of PDFs to zero in those regions.
176///
177/// After each PDF multiplication, the running product is compared with
178/// the cutOff parameter. If the running product is smaller than the
179/// cutOff value, the product series is terminated and remaining PDFs
180/// are not evaluated.
181///
182/// There is no magic value of the cutOff, the user should experiment
183/// to find the appropriate balance between speed and precision.
184/// If a cutoff is specified, the PDFs most likely to be small should
185/// be put first in the product. The default cutOff value is zero.
186
187RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgList& inPdfList, double cutOff) :
188 RooAbsPdf(name,title),
189 _cacheMgr(this,10),
190 _genCode(10),
191 _cutOff(cutOff),
192 _pdfList("!pdfs","List of PDFs",this),
193 _extendedIndex(-1),
194 _useDefaultGen(false),
195 _refRangeName(0),
196 _selfNorm(true)
197{
198 Int_t numExtended(0) ;
199 for(RooAbsArg * arg : inPdfList) {
200 RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
201 if (!pdf) {
202 coutW(InputArguments) << "RooProdPdf::RooProdPdf(" << GetName() << ") list arg "
203 << arg->GetName() << " is not a PDF, ignored" << endl ;
204 continue ;
205 }
206 _pdfList.add(*pdf) ;
207
208 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>("nset")) ;
209
210 if (pdf->canBeExtended()) {
212 numExtended++ ;
213 }
214 }
215
216 // Protect against multiple extended terms
217 if (numExtended>1) {
218 coutW(InputArguments) << "RooProdPdf::RooProdPdf(" << GetName()
219 << ") WARNING: multiple components with extended terms detected,"
220 << " product will not be extendible." << endl ;
221 _extendedIndex = -1 ;
222 }
223
225}
226
227
228
229////////////////////////////////////////////////////////////////////////////////
230/// Constructor from named argument list.
231/// \param[in] name Name used by RooFit
232/// \param[in] title Title used for plotting
233/// \param[in] fullPdfSet Set of "regular" PDFs that are normalised over all their observables
234/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 Optional arguments according to table below.
235///
236/// <table>
237/// <tr><th> Argument <th> Description
238/// <tr><td> `Conditional(pdfSet,depSet,depsAreCond=false)` <td> Add PDF to product with condition that it
239/// only be normalized over specified observables. Any remaining observables will be conditional observables.
240/// (Setting `depsAreCond` to true inverts this, so the observables in depSet will be the conditional observables.)
241/// </table>
242///
243/// For example, given a PDF \f$ F(x,y) \f$ and \f$ G(y) \f$,
244///
245/// `RooProdPdf("P", "P", G, Conditional(F,x))` will construct a 2-dimensional PDF as follows:
246/// \f[
247/// P(x,y) = \frac{G(y)}{\int_y G(y)} \cdot \frac{F(x,y)}{\int_x F(x,y)},
248/// \f]
249///
250/// which is a well normalised and properly defined PDF, but different from
251/// \f[
252/// P'(x,y) = \frac{F(x,y) \cdot G(y)}{\int_x\int_y F(x,y) \cdot G(y)}.
253/// \f]
254///
255/// In the former case, the \f$ y \f$ distribution of \f$ P \f$ is identical to that of \f$ G \f$, while
256/// \f$ F \f$ only is used to determine the correlation between \f$ X \f$ and \f$ Y \f$. In the latter
257/// case, the \f$ Y \f$ distribution is defined by the product of \f$ F \f$ and \f$ G \f$.
258///
259/// This \f$ P(x,y) \f$ construction is analoguous to generating events from \f$ F(x,y) \f$ with
260/// a prototype dataset sampled from \f$ G(y) \f$.
261
262RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet,
263 const RooCmdArg& arg1, const RooCmdArg& arg2,
264 const RooCmdArg& arg3, const RooCmdArg& arg4,
265 const RooCmdArg& arg5, const RooCmdArg& arg6,
266 const RooCmdArg& arg7, const RooCmdArg& arg8) :
267 RooAbsPdf(name,title),
268 _cacheMgr(this,10),
269 _genCode(10),
270 _cutOff(0),
271 _pdfList("!pdfs","List of PDFs",this),
272 _extendedIndex(-1),
273 _useDefaultGen(false),
274 _refRangeName(0),
275 _selfNorm(true)
276{
278 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
279 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
280 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
281 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
282
283 initializeFromCmdArgList(fullPdfSet,l) ;
285}
286
287
288
289////////////////////////////////////////////////////////////////////////////////
290/// Constructor from named argument list
291
292RooProdPdf::RooProdPdf(const char* name, const char* title,
293 const RooCmdArg& arg1, const RooCmdArg& arg2,
294 const RooCmdArg& arg3, const RooCmdArg& arg4,
295 const RooCmdArg& arg5, const RooCmdArg& arg6,
296 const RooCmdArg& arg7, const RooCmdArg& arg8) :
297 RooAbsPdf(name,title),
298 _cacheMgr(this,10),
299 _genCode(10),
300 _cutOff(0),
301 _pdfList("!pdfList","List of PDFs",this),
302 _extendedIndex(-1),
303 _useDefaultGen(false),
304 _refRangeName(0),
305 _selfNorm(true)
306{
308 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
309 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
310 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
311 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
312
315}
316
317
318
319////////////////////////////////////////////////////////////////////////////////
320/// Internal constructor from list of named arguments
321
322RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet, const RooLinkedList& cmdArgList) :
323 RooAbsPdf(name,title),
324 _cacheMgr(this,10),
325 _genCode(10),
326 _cutOff(0),
327 _pdfList("!pdfs","List of PDFs",this),
328 _extendedIndex(-1),
329 _useDefaultGen(false),
330 _refRangeName(0),
331 _selfNorm(true)
332{
333 initializeFromCmdArgList(fullPdfSet, cmdArgList) ;
335}
336
337
338
339////////////////////////////////////////////////////////////////////////////////
340/// Copy constructor
341
342RooProdPdf::RooProdPdf(const RooProdPdf& other, const char* name) :
343 RooAbsPdf(other,name),
344 _cacheMgr(other._cacheMgr,this),
345 _genCode(other._genCode),
346 _cutOff(other._cutOff),
347 _pdfList("!pdfs",this,other._pdfList),
348 _extendedIndex(other._extendedIndex),
349 _useDefaultGen(other._useDefaultGen),
350 _refRangeName(other._refRangeName),
351 _selfNorm(other._selfNorm),
352 _defNormSet(other._defNormSet)
353{
354 // Clone contents of normalizarion set list
355 for(auto const& nset : other._pdfNSetList) {
356 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>(nset->GetName()));
357 nset->snapshot(*_pdfNSetList.back());
358 }
360}
361
362
363
364////////////////////////////////////////////////////////////////////////////////
365/// Initialize RooProdPdf configuration from given list of RooCmdArg configuration arguments
366/// and set of 'regular' p.d.f.s in product
367
369{
370 Int_t numExtended(0) ;
371
372 // Process set of full PDFS
373 for(auto const* pdf : static_range_cast<RooAbsPdf*>(fullPdfSet)) {
374 _pdfList.add(*pdf) ;
375 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>("nset")) ;
376
377 if (pdf->canBeExtended()) {
379 numExtended++ ;
380 }
381
382 }
383
384 // Process list of conditional PDFs
385 for(auto * carg : static_range_cast<RooCmdArg*>(l)) {
386
387 if (0 == strcmp(carg->GetName(), "Conditional")) {
388
389 Int_t argType = carg->getInt(0) ;
390 RooArgSet* pdfSet = (RooArgSet*) carg->getSet(0) ;
391 RooArgSet* normSet = (RooArgSet*) carg->getSet(1) ;
392
393 for(auto * thePdf : static_range_cast<RooAbsPdf*>(*pdfSet)) {
394 _pdfList.add(*thePdf) ;
395
396 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>(0 == argType ? "nset" : "cset"));
397 normSet->snapshot(*_pdfNSetList.back());
398
399 if (thePdf->canBeExtended()) {
400 _extendedIndex = _pdfList.index(thePdf) ;
401 numExtended++ ;
402 }
403
404 }
405
406 } else if (0 != strlen(carg->GetName())) {
407 coutW(InputArguments) << "Unknown arg: " << carg->GetName() << endl ;
408 }
409 }
410
411 // Protect against multiple extended terms
412 if (numExtended>1) {
413 coutW(InputArguments) << "RooProdPdf::RooProdPdf(" << GetName()
414 << ") WARNING: multiple components with extended terms detected,"
415 << " product will not be extendible." << endl ;
416 _extendedIndex = -1 ;
417 }
418
419
420}
421
422
423
424////////////////////////////////////////////////////////////////////////////////
425/// Destructor
426
428{
430}
431
432
434 int code ;
435 auto cache = static_cast<CacheElem*>(_cacheMgr.getObj(nset, 0, &code)) ;
436
437 // If cache doesn't have our configuration, recalculate here
438 if (!cache) {
439 code = getPartIntList(nset, nullptr) ;
440 cache = static_cast<CacheElem*>(_cacheMgr.getObj(nset, 0, &code)) ;
441 }
442 return cache;
443}
444
445
446////////////////////////////////////////////////////////////////////////////////
447/// Calculate current value of object
448
450{
452}
453
454
455
456////////////////////////////////////////////////////////////////////////////////
457/// Calculate running product of pdfs terms, using the supplied
458/// normalization set in 'normSetList' for each component
459
460double RooProdPdf::calculate(const RooProdPdf::CacheElem& cache, bool /*verbose*/) const
461{
462 if (cache._isRearranged) {
463 if (dologD(Eval)) {
464 cxcoutD(Eval) << "RooProdPdf::calculate(" << GetName() << ") rearranged product calculation"
465 << " calculate: num = " << cache._rearrangedNum->GetName() << " = " << cache._rearrangedNum->getVal() << endl ;
466// cache._rearrangedNum->printComponentTree("",0,5) ;
467 cxcoutD(Eval) << "calculate: den = " << cache._rearrangedDen->GetName() << " = " << cache._rearrangedDen->getVal() << endl ;
468// cache._rearrangedDen->printComponentTree("",0,5) ;
469 }
470
471 return cache._rearrangedNum->getVal() / cache._rearrangedDen->getVal();
472 } else {
473
474 double value = 1.0;
475 assert(cache._normList.size() == cache._partList.size());
476 for (std::size_t i = 0; i < cache._partList.size(); ++i) {
477 const auto& partInt = static_cast<const RooAbsReal&>(cache._partList[i]);
478 const auto normSet = cache._normList[i].get();
479
480 const double piVal = partInt.getVal(normSet->getSize() > 0 ? normSet : nullptr);
481 value *= piVal ;
482 if (value <= _cutOff) break;
483 }
484
485 return value ;
486 }
487}
488
489////////////////////////////////////////////////////////////////////////////////
490/// Evaluate product of PDFs in batch mode.
491void RooProdPdf::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
492{
494 for (const RooAbsArg* i:_pdfList) {
495 auto span = dataMap.at(i);
496 // If the pdf doesn't depend on any observable (detected by it getting evaluated in scalar mode),
497 // it corresponds to a parameter constraint and should not be evaluated.
498 // These pdfs are evaluated in the RooConstraintSum that gets added to the likelihood in the end.
499 if(span.size() == 1) continue;
500 pdfs.push_back(span);
501 }
502 RooBatchCompute::ArgVector special{ static_cast<double>(pdfs.size()) };
504 dispatch->compute(stream, RooBatchCompute::ProdPdf, output, nEvents, pdfs, special);
505}
506
507namespace {
508
509template<class T>
510void eraseNullptrs(std::vector<T*>& v) {
511 v.erase(std::remove_if(v.begin(), v.end(), [](T* x){ return x == nullptr; } ), v.end());
512}
513
514void removeCommon(std::vector<RooAbsArg*> &v, std::span<RooAbsArg * const> other) {
515
516 for (auto const& arg : other) {
517 auto namePtrMatch = [&arg](const RooAbsArg* elm) {
518 return elm != nullptr && elm->namePtr() == arg->namePtr();
519 };
520
521 auto found = std::find_if(v.begin(), v.end(), namePtrMatch);
522 if(found != v.end()) {
523 *found = nullptr;
524 }
525 }
526 eraseNullptrs(v);
527}
528
529void addCommon(std::vector<RooAbsArg*> &v, std::vector<RooAbsArg*> const& o1, std::vector<RooAbsArg*> const& o2) {
530
531 for (auto const& arg : o1) {
532 auto namePtrMatch = [&arg](const RooAbsArg* elm) {
533 return elm->namePtr() == arg->namePtr();
534 };
535
536 if(std::find_if(o2.begin(), o2.end(), namePtrMatch) != o2.end()) {
537 v.push_back(arg);
538 }
539 }
540}
541
542}
543
544
545////////////////////////////////////////////////////////////////////////////////
546/// Factorize product in irreducible terms for given choice of integration/normalization
547
548void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& intSet,
549 RooLinkedList& termList, RooLinkedList& normList,
550 RooLinkedList& impDepList, RooLinkedList& crossDepList,
551 RooLinkedList& intList) const
552{
553 // List of all term dependents: normalization and imported
554 std::vector<RooArgSet> depAllList;
555 std::vector<RooArgSet> depIntNoNormList;
556
557 // Setup lists for factorization terms and their dependents
558 RooArgSet* term(0);
559 RooArgSet* termNormDeps(0);
560 RooArgSet* termIntDeps(0);
561 RooArgSet* termIntNoNormDeps(0);
562
563 std::vector<RooAbsArg*> pdfIntNoNormDeps;
564 std::vector<RooAbsArg*> pdfIntSet;
565 std::vector<RooAbsArg*> pdfNSet;
566 std::vector<RooAbsArg*> pdfCSet;
567 std::vector<RooAbsArg*> pdfNormDeps; // Dependents to be normalized for the PDF
568 std::vector<RooAbsArg*> pdfAllDeps; // All dependents of this PDF
569
570 // Loop over the PDFs
571 for(std::size_t iPdf = 0; iPdf < _pdfList.size(); ++iPdf) {
572 RooAbsPdf& pdf = static_cast<RooAbsPdf&>(_pdfList[iPdf]);
573 RooArgSet& pdfNSetOrig = *_pdfNSetList[iPdf];
574
575 pdfNSet.clear();
576 pdfCSet.clear();
577
578 // Make iterator over tree leaf node list to get the observables.
579 // This code is borrowed from RooAgsPdf::getObservables.
580 // RooAbsArg::treeNodeServer list is relatively expensive, so we only do it
581 // once and use it in a lambda function.
582 RooArgSet pdfLeafList("leafNodeServerList") ;
583 pdf.treeNodeServerList(&pdfLeafList,0,false,true,true) ;
584 auto getObservablesOfCurrentPdf = [&pdfLeafList](
585 std::vector<RooAbsArg*> & out,
586 const RooArgSet& dataList) {
587 for (const auto arg : pdfLeafList) {
588 if (arg->dependsOnValue(dataList) && arg->isLValue()) {
589 out.push_back(arg) ;
590 }
591 }
592 };
593
594 // Reduce pdfNSet to actual dependents
595 if (0 == strcmp("nset", pdfNSetOrig.GetName())) {
596 getObservablesOfCurrentPdf(pdfNSet, pdfNSetOrig);
597 } else if (0 == strcmp("cset", pdfNSetOrig.GetName())) {
598 getObservablesOfCurrentPdf(pdfNSet, normSet);
599 removeCommon(pdfNSet, pdfNSetOrig.get());
600 pdfCSet = pdfNSetOrig.get();
601 } else {
602 // Legacy mode. Interpret at NSet for backward compatibility
603 getObservablesOfCurrentPdf(pdfNSet, pdfNSetOrig);
604 }
605
606
607 pdfNormDeps.clear();
608 pdfAllDeps.clear();
609
610 // Make list of all dependents of this PDF
611 getObservablesOfCurrentPdf(pdfAllDeps, normSet);
612
613
614// cout << GetName() << ": pdf = " << pdf->GetName() << " pdfAllDeps = " << pdfAllDeps << " pdfNSet = " << *pdfNSet << " pdfCSet = " << *pdfCSet << endl;
615
616 // Make list of normalization dependents for this PDF;
617 if (!pdfNSet.empty()) {
618 // PDF is conditional
619 addCommon(pdfNormDeps, pdfAllDeps, pdfNSet);
620 } else {
621 // PDF is regular
622 pdfNormDeps = pdfAllDeps;
623 }
624
625// cout << GetName() << ": pdfNormDeps for " << pdf->GetName() << " = " << pdfNormDeps << endl;
626
627 pdfIntSet.clear();
628 getObservablesOfCurrentPdf(pdfIntSet, intSet) ;
629
630 // WVE if we have no norm deps, conditional observables should be taken out of pdfIntSet
631 if (pdfNormDeps.empty() && !pdfCSet.empty()) {
632 removeCommon(pdfIntSet, pdfCSet);
633// cout << GetName() << ": have no norm deps, removing conditional observables from intset" << endl;
634 }
635
636 pdfIntNoNormDeps.clear();
637 pdfIntNoNormDeps = pdfIntSet;
638 removeCommon(pdfIntNoNormDeps, pdfNormDeps);
639
640// cout << GetName() << ": pdf = " << pdf->GetName() << " intset = " << *pdfIntSet << " pdfIntNoNormDeps = " << pdfIntNoNormDeps << endl;
641
642 // Check if this PDF has dependents overlapping with one of the existing terms
643 bool done = false;
644 int j = 0;
645 auto lIter = termList.begin();
646 auto ldIter = normList.begin();
647 for(;lIter != termList.end(); (++lIter, ++ldIter, ++j)) {
648 termNormDeps = static_cast<RooArgSet*>(*ldIter);
649 term = static_cast<RooArgSet*>(*lIter);
650 // PDF should be added to existing term if
651 // 1) It has overlapping normalization dependents with any other PDF in existing term
652 // 2) It has overlapping dependents of any class for which integration is requested
653 // 3) If normalization happens over multiple ranges, and those ranges are both defined
654 // in either observable
655
656 bool normOverlap = termNormDeps->overlaps(pdfNormDeps.begin(), pdfNormDeps.end());
657 //bool intOverlap = pdfIntSet->overlaps(*termAllDeps);
658
659 if (normOverlap) {
660// cout << GetName() << ": this term overlaps with term " << (*term) << " in normalization observables" << endl;
661
662 term->add(pdf);
663 termNormDeps->add(pdfNormDeps.begin(), pdfNormDeps.end(), false);
664 depAllList[j].add(pdfAllDeps.begin(), pdfAllDeps.end(), false);
665 if (termIntDeps) {
666 termIntDeps->add(pdfIntSet.begin(), pdfIntSet.end(), false);
667 }
668 if (termIntNoNormDeps) {
669 termIntNoNormDeps->add(pdfIntNoNormDeps.begin(), pdfIntNoNormDeps.end(), false);
670 }
671 termIntNoNormDeps->add(pdfIntNoNormDeps.begin(), pdfIntNoNormDeps.end(), false);
672 done = true;
673 }
674 }
675
676 // If not, create a new term
677 if (!done) {
678 if (!(pdfNormDeps.empty() && pdfAllDeps.empty() &&
679 pdfIntSet.empty()) || 0 == normSet.getSize()) {
680// cout << GetName() << ": creating new term" << endl;
681 term = new RooArgSet("term");
682 termNormDeps = new RooArgSet("termNormDeps");
683 depAllList.emplace_back(pdfAllDeps.begin(), pdfAllDeps.end(), "termAllDeps");
684 termIntDeps = new RooArgSet(pdfIntSet.begin(), pdfIntSet.end(), "termIntDeps");
685 depIntNoNormList.emplace_back(pdfIntNoNormDeps.begin(), pdfIntNoNormDeps.end(), "termIntNoNormDeps");
686 termIntNoNormDeps = &depIntNoNormList.back();
687
688 term->add(pdf);
689 termNormDeps->add(pdfNormDeps.begin(), pdfNormDeps.end(), false);
690
691 termList.Add(term);
692 normList.Add(termNormDeps);
693 intList.Add(termIntDeps);
694 }
695 }
696
697 }
698
699 // Loop over list of terms again to determine 'imported' observables
700 int i = 0;
701 RooArgSet *normDeps;
702 auto lIter = termList.begin();
703 auto ldIter = normList.begin();
704 for(;lIter != termList.end(); (++lIter, ++ldIter, ++i)) {
705 normDeps = static_cast<RooArgSet*>(*ldIter);
706 term = static_cast<RooArgSet*>(*lIter);
707 // Make list of wholly imported dependents
708 RooArgSet impDeps(depAllList[i]);
709 impDeps.remove(*normDeps, true, true);
710 impDepList.Add(impDeps.snapshot());
711// cout << GetName() << ": list of imported dependents for term " << (*term) << " set to " << impDeps << endl ;
712
713 // Make list of cross dependents (term is self contained for these dependents,
714 // but components import dependents from other components)
715 auto crossDeps = std::unique_ptr<RooAbsCollection>{depIntNoNormList[i].selectCommon(*normDeps)};
716 crossDepList.Add(crossDeps->snapshot());
717// cout << GetName() << ": list of cross dependents for term " << (*term) << " set to " << *crossDeps << endl ;
718 }
719
720 return;
721}
722
723
724
725
726////////////////////////////////////////////////////////////////////////////////
727/// Return list of (partial) integrals of product terms for integration
728/// of p.d.f over observables iset while normalization over observables nset.
729/// Also return list of normalization sets to be used to evaluate
730/// each component in the list correctly.
731
732Int_t RooProdPdf::getPartIntList(const RooArgSet* nset, const RooArgSet* iset, const char* isetRangeName) const
733{
734// cout << " FOLKERT::RooProdPdf::getPartIntList(" << GetName() <<") nset = " << (nset?*nset:RooArgSet()) << endl
735// << " _normRange = " << _normRange << endl
736// << " iset = " << (iset?*iset:RooArgSet()) << endl
737// << " isetRangeName = " << (isetRangeName?isetRangeName:"<null>") << endl ;
738
739 // Check if this configuration was created before
740 Int_t sterileIdx(-1);
741
742 CacheElem* cache = (CacheElem*) _cacheMgr.getObj(nset,iset,&sterileIdx,isetRangeName);
743 if (cache) {
744 return _cacheMgr.lastIndex();
745 }
746
747 // Create containers for partial integral components to be generated
748 cache = new CacheElem;
749
750 // Factorize the product in irreducible terms for this nset
751 RooLinkedList terms, norms, imp, ints, cross;
752// cout << "RooProdPdf::getPIL -- now calling factorizeProduct()" << endl ;
753
754
755 // Normalization set used for factorization
756 RooArgSet factNset(nset ? (*nset) : _defNormSet);
757// cout << GetName() << "factNset = " << factNset << endl ;
758
759 factorizeProduct(factNset, iset ? (*iset) : RooArgSet(), terms, norms, imp, cross, ints);
760
761 RooArgSet *norm, *integ, *xdeps, *imps;
762
763 // Group irriducible terms that need to be (partially) integrated together
764 std::list<std::vector<RooArgSet*>> groupedList;
765 RooArgSet outerIntDeps;
766// cout << "RooProdPdf::getPIL -- now calling groupProductTerms()" << endl;
767 groupProductTerms(groupedList, outerIntDeps, terms, norms, imp, ints, cross);
768
769 // Loop over groups
770// cout<<"FK: pdf("<<GetName()<<") Starting selecting F(x|y)!"<<endl;
771 // Find groups of type F(x|y), i.e. termImpSet!=0, construct ratio object
772 map<string, RooArgSet> ratioTerms;
773 for (auto const& group : groupedList) {
774 if (1 == group.size()) {
775// cout<<"FK: Starting Single Term"<<endl;
776
777 RooArgSet* term = group[0];
778
779 Int_t termIdx = terms.IndexOf(term);
780 norm=(RooArgSet*) norms.At(termIdx);
781 imps=(RooArgSet*)imp.At(termIdx);
782 RooArgSet termNSet(*norm), termImpSet(*imps);
783
784// cout<<"FK: termImpSet.getSize() = "<<termImpSet.getSize()<< " " << termImpSet << endl;
785// cout<<"FK: _refRangeName = "<<_refRangeName<<endl;
786
787 if (termImpSet.getSize() > 0 && 0 != _refRangeName) {
788
789// cout << "WVE now here" << endl;
790
791 // WVE we can skip this if the ref range is equal to the normalization range
792 bool rangeIdentical(true);
793 RooFIter niter = termNSet.fwdIterator();
794 RooRealVar* normObs;
795// cout << "_normRange = " << _normRange << " _refRangeName = " << RooNameReg::str(_refRangeName) << endl ;
796 while ((normObs = (RooRealVar*) niter.next())) {
797 //FK: Here the refRange should be compared to _normRange, if it's set, and to the normObs range if it's not set
798 if (_normRange.Length() > 0) {
799 if (normObs->getMin(_normRange.Data()) != normObs->getMin(RooNameReg::str(_refRangeName))) rangeIdentical = false;
800 if (normObs->getMax(_normRange.Data()) != normObs->getMax(RooNameReg::str(_refRangeName))) rangeIdentical = false;
801 }
802 else{
803 if (normObs->getMin() != normObs->getMin(RooNameReg::str(_refRangeName))) rangeIdentical = false;
804 if (normObs->getMax() != normObs->getMax(RooNameReg::str(_refRangeName))) rangeIdentical = false;
805 }
806 }
807// cout<<"FK: rangeIdentical Single = "<<(rangeIdentical ? 'T':'F')<<endl;
808 // coverity[CONSTANT_EXPRESSION_RESULT]
809 // LM : avoid making integral ratio if range is the same. Why was not included ??? (same at line 857)
810 if (!rangeIdentical ) {
811// cout << "PREPARING RATIO HERE (SINGLE TERM)" << endl ;
812 RooAbsReal* ratio = makeCondPdfRatioCorr(*(RooAbsReal*)term->first(), termNSet, termImpSet, normRange(), RooNameReg::str(_refRangeName));
813 ostringstream str; termImpSet.printValue(str);
814// cout << GetName() << "inserting ratio term" << endl;
815 ratioTerms[str.str()].add(*ratio);
816 }
817 }
818
819 } else {
820// cout<<"FK: Starting Composite Term"<<endl;
821
822 RooArgSet compTermSet, compTermNorm;
823 for (auto const& term : group) {
824
825 Int_t termIdx = terms.IndexOf(term);
826 norm=(RooArgSet*) norms.At(termIdx);
827 imps=(RooArgSet*)imp.At(termIdx);
828 RooArgSet termNSet(*norm), termImpSet(*imps);
829
830 if (termImpSet.getSize() > 0 && 0 != _refRangeName) {
831
832 // WVE we can skip this if the ref range is equal to the normalization range
833 bool rangeIdentical(true);
834 RooFIter niter = termNSet.fwdIterator();
835 RooRealVar* normObs;
836 //FK: Here the refRange should be compared to _normRange, if it's set, and to the normObs range if it's not set
837 if(_normRange.Length() > 0) {
838 while ((normObs = (RooRealVar*) niter.next())) {
839 if (normObs->getMin(_normRange.Data()) != normObs->getMin(RooNameReg::str(_refRangeName))) rangeIdentical = false;
840 if (normObs->getMax(_normRange.Data()) != normObs->getMax(RooNameReg::str(_refRangeName))) rangeIdentical = false;
841 }
842 } else {
843 while ((normObs = (RooRealVar*) niter.next())) {
844 if (normObs->getMin() != normObs->getMin(RooNameReg::str(_refRangeName))) rangeIdentical = false;
845 if (normObs->getMax() != normObs->getMax(RooNameReg::str(_refRangeName))) rangeIdentical = false;
846 }
847 }
848// cout<<"FK: rangeIdentical Composite = "<<(rangeIdentical ? 'T':'F') <<endl;
849 if (!rangeIdentical ) {
850// cout << "PREPARING RATIO HERE (COMPOSITE TERM)" << endl ;
851 RooAbsReal* ratio = makeCondPdfRatioCorr(*(RooAbsReal*)term->first(), termNSet, termImpSet, normRange(), RooNameReg::str(_refRangeName));
852 ostringstream str; termImpSet.printValue(str);
853 ratioTerms[str.str()].add(*ratio);
854 }
855 }
856 }
857 }
858
859 }
860
861 // Find groups with y as termNSet
862 // Replace G(y) with (G(y),ratio)
863 for (auto const& group : groupedList) {
864 if (1 == group.size()) {
865 RooArgSet* term = group[0];
866
867 Int_t termIdx = terms.IndexOf(term);
868 norm = (RooArgSet*) norms.At(termIdx);
869 imps = (RooArgSet*) imp.At(termIdx);
870 RooArgSet termNSet(*norm), termImpSet(*imps);
871
872 // If termNset matches index of ratioTerms, insert ratio here
873 ostringstream str; termNSet.printValue(str);
874 if (ratioTerms[str.str()].getSize() > 0) {
875// cout << "MUST INSERT RATIO OBJECT IN TERM (SINGLE) " << *term << endl;
876 term->add(ratioTerms[str.str()]);
877 }
878 } else {
879 RooArgSet compTermSet, compTermNorm;
880 for (auto const& term : group) {
881 Int_t termIdx = terms.IndexOf(term);
882 norm = (RooArgSet*) norms.At(termIdx);
883 imps = (RooArgSet*) imp.At(termIdx);
884 RooArgSet termNSet(*norm), termImpSet(*imps);
885
886 // If termNset matches index of ratioTerms, insert ratio here
887 ostringstream str; termNSet.printValue(str);
888 if (ratioTerms[str.str()].getSize() > 0) {
889// cout << "MUST INSERT RATIO OBJECT IN TERM (COMPOSITE)" << *term << endl;
890 term->add(ratioTerms[str.str()]);
891 }
892 }
893 }
894 }
895
896 for (auto const& group : groupedList) {
897// cout << GetName() << ":now processing group" << endl;
898// group->Print("1");
899
900 if (1 == group.size()) {
901// cout << "processing atomic item" << endl;
902 RooArgSet* term = group[0];
903
904 Int_t termIdx = terms.IndexOf(term);
905 norm = (RooArgSet*) norms.At(termIdx);
906 integ = (RooArgSet*) ints.At(termIdx);
907 xdeps = (RooArgSet*) cross.At(termIdx);
908 imps = (RooArgSet*) imp.At(termIdx);
909
910 RooArgSet termNSet, termISet, termXSet, termImpSet;
911
912 // Take list of normalization, integrated dependents from factorization algorithm
913 termISet.add(*integ);
914 termNSet.add(*norm);
915
916 // Cross-imported integrated dependents
917 termXSet.add(*xdeps);
918 termImpSet.add(*imps);
919
920// cout << GetName() << ": termISet = " << termISet << endl;
921// cout << GetName() << ": termNSet = " << termNSet << endl;
922// cout << GetName() << ": termXSet = " << termXSet << endl;
923// cout << GetName() << ": termImpSet = " << termImpSet << endl;
924
925 // Add prefab term to partIntList.
926 bool isOwned(false);
927 vector<RooAbsReal*> func = processProductTerm(nset, iset, isetRangeName, term, termNSet, termISet, isOwned);
928 if (func[0]) {
929 cache->_partList.add(*func[0]);
930 if (isOwned) cache->_ownedList.addOwned(*func[0]);
931
932 cache->_normList.emplace_back(norm->snapshot(false));
933
934 cache->_numList.addOwned(*func[1]);
935 cache->_denList.addOwned(*func[2]);
936// cout << "func[0]=" << func[0]->IsA()->GetName() << "::" << func[0]->GetName() << endl;
937// cout << "func[1]=" << func[1]->IsA()->GetName() << "::" << func[1]->GetName() << endl;
938// cout << "func[2]=" << func[2]->IsA()->GetName() << "::" << func[2]->GetName() << endl;
939 }
940 } else {
941// cout << "processing composite item" << endl;
942 RooArgSet compTermSet, compTermNorm, compTermNum, compTermDen;
943 for (auto const& term : group) {
944// cout << GetName() << ": processing term " << (*term) << " of composite item" << endl ;
945 Int_t termIdx = terms.IndexOf(term);
946 norm = (RooArgSet*) norms.At(termIdx);
947 integ = (RooArgSet*) ints.At(termIdx);
948 xdeps = (RooArgSet*) cross.At(termIdx);
949 imps = (RooArgSet*) imp.At(termIdx);
950
951 RooArgSet termNSet, termISet, termXSet, termImpSet;
952 termISet.add(*integ);
953 termNSet.add(*norm);
954 termXSet.add(*xdeps);
955 termImpSet.add(*imps);
956
957 // Remove outer integration dependents from termISet
958 termISet.remove(outerIntDeps, true, true);
959// cout << "termISet = "; termISet.Print("1");
960
961// cout << GetName() << ": termISet = " << termISet << endl;
962// cout << GetName() << ": termNSet = " << termNSet << endl;
963// cout << GetName() << ": termXSet = " << termXSet << endl;
964// cout << GetName() << ": termImpSet = " << termImpSet << endl;
965
966 bool isOwned = false;
967 vector<RooAbsReal*> func = processProductTerm(nset, iset, isetRangeName, term, termNSet, termISet, isOwned, true);
968// cout << GetName() << ": created composite term component " << func[0]->GetName() << endl;
969 if (func[0]) {
970 compTermSet.add(*func[0]);
971 if (isOwned) cache->_ownedList.addOwned(*func[0]);
972 compTermNorm.add(*norm, false);
973
974 compTermNum.add(*func[1]);
975 compTermDen.add(*func[2]);
976 //cache->_numList.add(*func[1]);
977 //cache->_denList.add(*func[2]);
978
979// cout << "func[0]=" << func[0]->IsA()->GetName() << "::" << func[0]->GetName() << endl;
980// cout << "func[1]=" << func[1]->IsA()->GetName() << "::" << func[1]->GetName() << endl;
981// cout << "func[2]=" << func[2]->IsA()->GetName() << "::" << func[2]->GetName() << endl;
982 }
983 }
984
985// cout << GetName() << ": constructing special composite product" << endl;
986// cout << GetName() << ": compTermSet = " ; compTermSet.Print("1");
987
988 // WVE THIS NEEDS TO BE REARRANGED
989
990 // compTermset is set van partial integrals to be multiplied
991 // prodtmp = product (compTermSet)
992 // inttmp = int ( prodtmp ) d (outerIntDeps) _range_isetRangeName
993
994 const std::string prodname = makeRGPPName("SPECPROD", compTermSet, outerIntDeps, RooArgSet(), isetRangeName);
995 RooProduct* prodtmp = new RooProduct(prodname.c_str(), prodname.c_str(), compTermSet);
996 cache->_ownedList.addOwned(*prodtmp);
997
998 const std::string intname = makeRGPPName("SPECINT", compTermSet, outerIntDeps, RooArgSet(), isetRangeName);
999 RooRealIntegral* inttmp = new RooRealIntegral(intname.c_str(), intname.c_str(), *prodtmp, outerIntDeps, 0, 0, isetRangeName);
1000 inttmp->setStringAttribute("PROD_TERM_TYPE", "SPECINT");
1001
1002 cache->_ownedList.addOwned(*inttmp);
1003 cache->_partList.add(*inttmp);
1004
1005 // Product of numerator terms
1006 const string prodname_num = makeRGPPName("SPECPROD_NUM", compTermNum, RooArgSet(), RooArgSet(), 0);
1007 RooProduct* prodtmp_num = new RooProduct(prodname_num.c_str(), prodname_num.c_str(), compTermNum);
1008 prodtmp_num->addOwnedComponents(compTermNum);
1009 cache->_ownedList.addOwned(*prodtmp_num);
1010
1011 // Product of denominator terms
1012 const string prodname_den = makeRGPPName("SPECPROD_DEN", compTermDen, RooArgSet(), RooArgSet(), 0);
1013 RooProduct* prodtmp_den = new RooProduct(prodname_den.c_str(), prodname_den.c_str(), compTermDen);
1014 prodtmp_den->addOwnedComponents(compTermDen);
1015 cache->_ownedList.addOwned(*prodtmp_den);
1016
1017 // Ratio
1018 string name = Form("SPEC_RATIO(%s,%s)", prodname_num.c_str(), prodname_den.c_str());
1019 RooFormulaVar* ndr = new RooFormulaVar(name.c_str(), "@0/@1", RooArgList(*prodtmp_num, *prodtmp_den));
1020
1021 // Integral of ratio
1022 RooAbsReal* numtmp = ndr->createIntegral(outerIntDeps,isetRangeName);
1023 numtmp->addOwnedComponents(*ndr);
1024
1025 cache->_numList.addOwned(*numtmp);
1027 cache->_normList.emplace_back(compTermNorm.snapshot(false));
1028 }
1029 }
1030
1031 // Store the partial integral list and return the assigned code
1032 Int_t returnCode = _cacheMgr.setObj(nset, iset, (RooAbsCacheElement*)cache, RooNameReg::ptr(isetRangeName));
1033
1034 // WVE DEBUG PRINTING
1035// cout << "RooProdPdf::getPartIntList(" << GetName() << ") made cache " << cache << " with the following nset pointers ";
1036// TIterator* nliter = nsetList->MakeIterator();
1037// RooArgSet* ns;
1038// while((ns=(RooArgSet*)nliter->Next())) {
1039// cout << ns << " ";
1040// }
1041// cout << endl;
1042// delete nliter;
1043
1044// cout << " FOLKERT::RooProdPdf::getPartIntList END(" << GetName() <<") nset = " << (nset?*nset:RooArgSet()) << endl
1045// << " _normRange = " << _normRange << endl
1046// << " iset = " << (iset?*iset:RooArgSet()) << endl
1047// << " partList = ";
1048// if(partListPointer) partListPointer->Print();
1049// cout << " nsetList = ";
1050// if(nsetListPointer) nsetListPointer->Print("");
1051// cout << " code = " << returnCode << endl
1052// << " isetRangeName = " << (isetRangeName?isetRangeName:"<null>") << endl;
1053
1054
1055 // Need to rearrange product in case of multiple ranges
1056 if (_normRange.Contains(",")) {
1057 rearrangeProduct(*cache);
1058 }
1059
1060 // We own contents of all lists filled by factorizeProduct()
1061 terms.Delete();
1062 ints.Delete();
1063 imp.Delete();
1064 norms.Delete();
1065 cross.Delete();
1066
1067 return returnCode;
1068}
1069
1070
1071
1072
1073////////////////////////////////////////////////////////////////////////////////
1074/// For single normalization ranges
1075
1076RooAbsReal* RooProdPdf::makeCondPdfRatioCorr(RooAbsReal& pdf, const RooArgSet& termNset, const RooArgSet& /*termImpSet*/, const char* normRangeTmp, const char* refRange) const
1077{
1078 RooAbsReal* ratio_num = pdf.createIntegral(termNset,normRangeTmp) ;
1079 RooAbsReal* ratio_den = pdf.createIntegral(termNset,refRange) ;
1080 RooFormulaVar* ratio = new RooFormulaVar(Form("ratio(%s,%s)",ratio_num->GetName(),ratio_den->GetName()),"@0/@1",
1081 RooArgList(*ratio_num,*ratio_den)) ;
1082
1083 ratio->addOwnedComponents(RooArgSet(*ratio_num,*ratio_den)) ;
1084 ratio->setAttribute("RATIO_TERM") ;
1085 return ratio ;
1086}
1087
1088
1089
1090
1091////////////////////////////////////////////////////////////////////////////////
1092
1094{
1095 RooAbsReal* part, *num, *den ;
1096 RooArgSet nomList ;
1097
1098 list<string> rangeComps ;
1099 {
1100 char* buf = new char[strlen(_normRange.Data()) + 1] ;
1101 strcpy(buf,_normRange.Data()) ;
1102 char* save(0) ;
1103 char* token = R__STRTOK_R(buf,",",&save) ;
1104 while(token) {
1105 rangeComps.push_back(token) ;
1106 token = R__STRTOK_R(0,",",&save) ;
1107 }
1108 delete[] buf;
1109 }
1110
1111
1112 map<string,RooArgSet> denListList ;
1113 RooArgSet specIntDeps ;
1114 string specIntRange ;
1115
1116// cout << "THIS IS REARRANGEPRODUCT" << endl ;
1117
1118 RooFIter iterp = cache._partList.fwdIterator() ;
1119 RooFIter iter1 = cache._numList.fwdIterator() ;
1120 RooFIter iter2 = cache._denList.fwdIterator() ;
1121 while((part=(RooAbsReal*)iterp.next())) {
1122
1123 num = (RooAbsReal*) iter1.next() ;
1124 den = (RooAbsReal*) iter2.next() ;
1125
1126// cout << "now processing part " << part->GetName() << " of type " << part->getStringAttribute("PROD_TERM_TYPE") << endl ;
1127// cout << "corresponding numerator = " << num->GetName() << endl ;
1128// cout << "corresponding denominator = " << den->GetName() << endl ;
1129
1130
1131 RooFormulaVar* ratio(0) ;
1132 RooArgSet origNumTerm ;
1133
1134 if (string("SPECINT")==part->getStringAttribute("PROD_TERM_TYPE")) {
1135
1136 RooRealIntegral* orig = (RooRealIntegral*) num;
1137 RooFormulaVar* specratio = (RooFormulaVar*) &orig->integrand() ;
1138 RooProduct* func = (RooProduct*) specratio->getParameter(0) ;
1139
1140 std::unique_ptr<RooArgSet> components{orig->getComponents()};
1141 for(RooAbsArg * carg : *components) {
1142 if (carg->getAttribute("RATIO_TERM")) {
1143 ratio = (RooFormulaVar*)carg ;
1144 break ;
1145 }
1146 }
1147
1148 if (ratio) {
1149 RooCustomizer cust(*func,"blah") ;
1150 cust.replaceArg(*ratio,RooFit::RooConst(1)) ;
1151 RooAbsArg* funcCust = cust.build() ;
1152// cout << "customized function = " << endl ;
1153// funcCust->printComponentTree() ;
1154 nomList.add(*funcCust) ;
1155 } else {
1156 nomList.add(*func) ;
1157 }
1158
1159
1160 } else {
1161
1162 // Find the ratio term
1163 RooAbsReal* func = num;
1164 // If top level object is integral, navigate to integrand
1165 if (func->InheritsFrom(RooRealIntegral::Class())) {
1166 func = (RooAbsReal*) &((RooRealIntegral*)(func))->integrand();
1167 }
1168 if (func->InheritsFrom(RooProduct::Class())) {
1169// cout << "product term found: " ; func->Print() ;
1170 for(RooAbsArg * arg : static_cast<RooProduct*>(func)->components()) {
1171 if (arg->getAttribute("RATIO_TERM")) {
1172 ratio = (RooFormulaVar*)(arg) ;
1173 } else {
1174 origNumTerm.add(*arg) ;
1175 }
1176 }
1177 }
1178
1179 if (ratio) {
1180// cout << "Found ratio term in numerator: " << ratio->GetName() << endl ;
1181// cout << "Adding only original term to numerator: " << origNumTerm << endl ;
1182 nomList.add(origNumTerm) ;
1183 } else {
1184 nomList.add(*num) ;
1185 }
1186
1187 }
1188
1189 for (list<string>::iterator iter = rangeComps.begin() ; iter != rangeComps.end() ; ++iter) {
1190 // If denominator is an integral, make a clone with the integration range adjusted to
1191 // the selected component of the normalization integral
1192// cout << "NOW PROCESSING DENOMINATOR " << den->IsA()->GetName() << "::" << den->GetName() << endl ;
1193
1194 if (string("SPECINT")==part->getStringAttribute("PROD_TERM_TYPE")) {
1195
1196// cout << "create integral: SPECINT case" << endl ;
1197 RooRealIntegral* orig = (RooRealIntegral*) num;
1198 RooFormulaVar* specRatio = (RooFormulaVar*) &orig->integrand() ;
1199 specIntDeps.add(orig->intVars()) ;
1200 if (orig->intRange()) {
1201 specIntRange = orig->intRange() ;
1202 }
1203 //RooProduct* numtmp = (RooProduct*) specRatio->getParameter(0) ;
1204 RooProduct* dentmp = (RooProduct*) specRatio->getParameter(1) ;
1205
1206// cout << "numtmp = " << numtmp->IsA()->GetName() << "::" << numtmp->GetName() << endl ;
1207// cout << "dentmp = " << dentmp->IsA()->GetName() << "::" << dentmp->GetName() << endl ;
1208
1209// cout << "denominator components are " << dentmp->components() << endl ;
1210 RooArgSet comps(dentmp->components()) ;
1211 RooFIter piter = comps.fwdIterator() ;
1212 RooAbsReal* parg ;
1213 while((parg=(RooAbsReal*)piter.next())) {
1214// cout << "now processing denominator component " << parg->IsA()->GetName() << "::" << parg->GetName() << endl ;
1215
1216 if (ratio && parg->dependsOn(*ratio)) {
1217// cout << "depends in value of ratio" << endl ;
1218
1219 // Make specialize ratio instance
1220 RooAbsReal* specializedRatio = specializeRatio(*(RooFormulaVar*)ratio,iter->c_str()) ;
1221
1222// cout << "specRatio = " << endl ;
1223// specializedRatio->printComponentTree() ;
1224
1225 // Replace generic ratio with specialized ratio
1226 RooAbsArg *partCust(0) ;
1227 if (parg->InheritsFrom(RooAddition::Class())) {
1228
1229
1230
1231 RooAddition* tmpadd = (RooAddition*)(parg) ;
1232
1233 RooCustomizer cust(*tmpadd->list1().first(),Form("blah_%s",iter->c_str())) ;
1234 cust.replaceArg(*ratio,*specializedRatio) ;
1235 partCust = cust.build() ;
1236
1237 } else {
1238 RooCustomizer cust(*parg,Form("blah_%s",iter->c_str())) ;
1239 cust.replaceArg(*ratio,*specializedRatio) ;
1240 partCust = cust.build() ;
1241 }
1242
1243 // Print customized denominator
1244// cout << "customized function = " << endl ;
1245// partCust->printComponentTree() ;
1246
1247 RooAbsReal* specializedPartCust = specializeIntegral(*(RooAbsReal*)partCust,iter->c_str()) ;
1248
1249 // Finally divide again by ratio
1250 string name = Form("%s_divided_by_ratio",specializedPartCust->GetName()) ;
1251 RooFormulaVar* specIntFinal = new RooFormulaVar(name.c_str(),"@0/@1",RooArgList(*specializedPartCust,*specializedRatio)) ;
1252
1253 denListList[*iter].add(*specIntFinal) ;
1254 } else {
1255
1256// cout << "does NOT depend on value of ratio" << endl ;
1257// parg->Print("t") ;
1258
1259 denListList[*iter].add(*specializeIntegral(*parg,iter->c_str())) ;
1260
1261 }
1262 }
1263// cout << "end iteration over denominator components" << endl ;
1264 } else {
1265
1266 if (ratio) {
1267
1268 RooAbsReal* specRatio = specializeRatio(*(RooFormulaVar*)ratio,iter->c_str()) ;
1269
1270 // If integral is 'Int r(y)*g(y) dy ' then divide a posteriori by r(y)
1271// cout << "have ratio, orig den = " << den->GetName() << endl ;
1272
1273 RooArgSet tmp(origNumTerm) ;
1274 tmp.add(*specRatio) ;
1275 const string pname = makeRGPPName("PROD",tmp,RooArgSet(),RooArgSet(),0) ;
1276 RooProduct* specDenProd = new RooProduct(pname.c_str(),pname.c_str(),tmp) ;
1277 RooAbsReal* specInt(0) ;
1278
1280 specInt = specDenProd->createIntegral(((RooRealIntegral*)den)->intVars(),iter->c_str()) ;
1281 } else if (den->InheritsFrom(RooAddition::Class())) {
1282 RooAddition* orig = (RooAddition*)den ;
1283 RooRealIntegral* origInt = (RooRealIntegral*) orig->list1().first() ;
1284 specInt = specDenProd->createIntegral(origInt->intVars(),iter->c_str()) ;
1285 } else {
1286 throw string("this should not happen") ;
1287 }
1288
1289 //RooAbsReal* specInt = specializeIntegral(*den,iter->c_str()) ;
1290 string name = Form("%s_divided_by_ratio",specInt->GetName()) ;
1291 RooFormulaVar* specIntFinal = new RooFormulaVar(name.c_str(),"@0/@1",RooArgList(*specInt,*specRatio)) ;
1292 denListList[*iter].add(*specIntFinal) ;
1293 } else {
1294 denListList[*iter].add(*specializeIntegral(*den,iter->c_str())) ;
1295 }
1296
1297 }
1298 }
1299
1300 }
1301
1302 // Do not rearrage terms if numerator and denominator are effectively empty
1303 if (nomList.getSize()==0) {
1304 return ;
1305 }
1306
1307 string name = Form("%s_numerator",GetName()) ;
1308 // WVE FIX THIS (2)
1309
1310 RooAbsReal* numerator = new RooProduct(name.c_str(),name.c_str(),nomList) ;
1311
1312 RooArgSet products ;
1313// cout << "nomList = " << nomList << endl ;
1314 for (map<string,RooArgSet>::iterator iter = denListList.begin() ; iter != denListList.end() ; ++iter) {
1315// cout << "denList[" << iter->first << "] = " << iter->second << endl ;
1316 name = Form("%s_denominator_comp_%s",GetName(),iter->first.c_str()) ;
1317 // WVE FIX THIS (2)
1318 RooProduct* prod_comp = new RooProduct(name.c_str(),name.c_str(),iter->second) ;
1319 products.add(*prod_comp) ;
1320 }
1321 name = Form("%s_denominator_sum",GetName()) ;
1322 RooAbsReal* norm = new RooAddition(name.c_str(),name.c_str(),products) ;
1323 norm->addOwnedComponents(products) ;
1324
1325 if (specIntDeps.getSize()>0) {
1326 // Apply posterior integration required for SPECINT case
1327
1328 string namesr = Form("SPEC_RATIO(%s,%s)",numerator->GetName(),norm->GetName()) ;
1329 RooFormulaVar* ndr = new RooFormulaVar(namesr.c_str(),"@0/@1",RooArgList(*numerator,*norm)) ;
1330
1331 // Integral of ratio
1332 RooAbsReal* numtmp = ndr->createIntegral(specIntDeps,specIntRange.c_str()) ;
1333
1334 numerator = numtmp ;
1335 norm = (RooAbsReal*) RooFit::RooConst(1).Clone() ;
1336 }
1337
1338
1339// cout << "numerator" << endl ;
1340// numerator->printComponentTree("",0,5) ;
1341// cout << "denominator" << endl ;
1342// norm->printComponentTree("",0,5) ;
1343
1344
1345 // WVE DEBUG
1346 //RooMsgService::instance().debugWorkspace()->import(RooArgSet(*numerator,*norm)) ;
1347
1348 cache._rearrangedNum.reset(numerator);
1349 cache._rearrangedDen.reset(norm);
1350 cache._isRearranged = true ;
1351
1352}
1353
1354
1355////////////////////////////////////////////////////////////////////////////////
1356
1357RooAbsReal* RooProdPdf::specializeRatio(RooFormulaVar& input, const char* targetRangeName) const
1358{
1359 RooRealIntegral* numint = (RooRealIntegral*) input.getParameter(0) ;
1360 RooRealIntegral* denint = (RooRealIntegral*) input.getParameter(1) ;
1361
1362 RooAbsReal* numint_spec = specializeIntegral(*numint,targetRangeName) ;
1363
1364 RooAbsReal* ret = new RooFormulaVar(Form("ratio(%s,%s)",numint_spec->GetName(),denint->GetName()),"@0/@1",RooArgList(*numint_spec,*denint)) ;
1365 ret->addOwnedComponents(*numint_spec) ;
1366
1367 return ret ;
1368}
1369
1370
1371
1372////////////////////////////////////////////////////////////////////////////////
1373
1374RooAbsReal* RooProdPdf::specializeIntegral(RooAbsReal& input, const char* targetRangeName) const
1375{
1376 if (input.InheritsFrom(RooRealIntegral::Class())) {
1377
1378 // If input is integral, recreate integral but override integration range to be targetRangeName
1380// cout << "creating integral: integrand = " << orig->integrand().GetName() << " vars = " << orig->intVars() << " range = " << targetRangeName << endl ;
1381 return orig->integrand().createIntegral(orig->intVars(),targetRangeName) ;
1382
1383 } else if (input.InheritsFrom(RooAddition::Class())) {
1384
1385 // If input is sum of integrals, recreate integral from first component of set, but override integration range to be targetRangeName
1386 RooAddition* orig = (RooAddition*)&input ;
1387 RooRealIntegral* origInt = (RooRealIntegral*) orig->list1().first() ;
1388// cout << "creating integral from addition: integrand = " << origInt->integrand().GetName() << " vars = " << origInt->intVars() << " range = " << targetRangeName << endl ;
1389 return origInt->integrand().createIntegral(origInt->intVars(),targetRangeName) ;
1390
1391 } else {
1392
1393// cout << "specializeIntegral: unknown input type " << input.IsA()->GetName() << "::" << input.GetName() << endl ;
1394 }
1395
1396 return &input ;
1397}
1398
1399
1400////////////////////////////////////////////////////////////////////////////////
1401/// Group product into terms that can be calculated independently
1402
1403void RooProdPdf::groupProductTerms(std::list<std::vector<RooArgSet*>>& groupedTerms, RooArgSet& outerIntDeps,
1404 const RooLinkedList& terms, const RooLinkedList& norms,
1405 const RooLinkedList& imps, const RooLinkedList& ints, const RooLinkedList& /*cross*/) const
1406{
1407 // Start out with each term in its own group
1408 for(auto * term : static_range_cast<RooArgSet*>(terms)) {
1409 groupedTerms.emplace_back();
1410 groupedTerms.back().emplace_back(term) ;
1411 }
1412
1413 // Make list of imported dependents that occur in any term
1414 RooArgSet allImpDeps ;
1415 for(auto * impDeps : static_range_cast<RooArgSet*>(imps)) {
1416 allImpDeps.add(*impDeps,false) ;
1417 }
1418
1419 // Make list of integrated dependents that occur in any term
1420 RooArgSet allIntDeps ;
1421 for(auto * intDeps : static_range_cast<RooArgSet*>(ints)) {
1422 allIntDeps.add(*intDeps,false) ;
1423 }
1424
1425 outerIntDeps.removeAll() ;
1426 outerIntDeps.add(*std::unique_ptr<RooArgSet>{static_cast<RooArgSet*>(allIntDeps.selectCommon(allImpDeps))});
1427
1428 // Now iteratively merge groups that should be (partially) integrated together
1429 for(RooAbsArg * outerIntDep : outerIntDeps) {
1430
1431 // Collect groups that feature this dependent
1432 std::vector<RooArgSet*>* newGroup = nullptr ;
1433
1434 // Loop over groups
1435 bool needMerge = false ;
1436 auto group = groupedTerms.begin();
1437 auto nGroups = groupedTerms.size();
1438 for (size_t iGroup = 0; iGroup < nGroups; ++iGroup) {
1439
1440 // See if any term in this group depends in any ay on outerDepInt
1441 for (auto const& term2 : *group) {
1442
1443 Int_t termIdx = terms.IndexOf(term2) ;
1444 RooArgSet* termNormDeps = (RooArgSet*) norms.At(termIdx) ;
1445 RooArgSet* termIntDeps = (RooArgSet*) ints.At(termIdx) ;
1446 RooArgSet* termImpDeps = (RooArgSet*) imps.At(termIdx) ;
1447
1448 if (termNormDeps->contains(*outerIntDep) ||
1449 termIntDeps->contains(*outerIntDep) ||
1450 termImpDeps->contains(*outerIntDep)) {
1451 needMerge = true ;
1452 }
1453
1454 }
1455
1456 if (needMerge) {
1457 // Create composite group if not yet existing
1458 if (newGroup==nullptr) {
1459 groupedTerms.emplace_back() ;
1460 newGroup = &groupedTerms.back() ;
1461 }
1462
1463 // Add terms of this group to new term
1464 for (auto& term2 : *group) {
1465 newGroup->emplace_back(term2) ;
1466 }
1467
1468 // Remove this group from list and delete it (but not its contents)
1469 group = groupedTerms.erase(group);
1470 } else {
1471 ++group;
1472 }
1473 }
1474
1475 }
1476}
1477
1478
1479
1480////////////////////////////////////////////////////////////////////////////////
1481/// Calculate integrals of factorized product terms over observables iset while normalized
1482/// to observables in nset.
1483
1484std::vector<RooAbsReal*> RooProdPdf::processProductTerm(const RooArgSet* nset, const RooArgSet* iset, const char* isetRangeName,
1485 const RooArgSet* term,const RooArgSet& termNSet, const RooArgSet& termISet,
1486 bool& isOwned, bool forceWrap) const
1487{
1488// cout << " FOLKERT::RooProdPdf(" << GetName() <<") processProductTerm nset = " << (nset?*nset:RooArgSet()) << endl
1489// << " _normRange = " << _normRange << endl
1490// << " iset = " << (iset?*iset:RooArgSet()) << endl
1491// << " isetRangeName = " << (isetRangeName?isetRangeName:"<null>") << endl
1492// << " term = " << (term?*term:RooArgSet()) << endl
1493// << " termNSet = " << termNSet << endl
1494// << " termISet = " << termISet << endl
1495// << " isOwned = " << isOwned << endl
1496// << " forceWrap = " << forceWrap << endl ;
1497
1498 vector<RooAbsReal*> ret(3) ; ret[0] = 0 ; ret[1] = 0 ; ret[2] = 0 ;
1499
1500 // CASE I: factorizing term: term is integrated over all normalizing observables
1501 // -----------------------------------------------------------------------------
1502 // Check if all observbales of this term are integrated. If so the term cancels
1503 if (termNSet.getSize()>0 && termNSet.getSize()==termISet.getSize() && isetRangeName==0) {
1504
1505
1506 //cout << "processProductTerm(" << GetName() << ") case I " << endl ;
1507
1508 // Term factorizes
1509 return ret ;
1510 }
1511
1512 // CASE II: Dropped terms: if term is entirely unnormalized, it should be dropped
1513 // ------------------------------------------------------------------------------
1514 if (nset && termNSet.getSize()==0) {
1515
1516 //cout << "processProductTerm(" << GetName() << ") case II " << endl ;
1517
1518 // Drop terms that are not asked to be normalized
1519 return ret ;
1520 }
1521
1522 if (iset && termISet.getSize()>0) {
1523 if (term->getSize()==1) {
1524
1525 // CASE IIIa: Normalized and partially integrated single PDF term
1526 //---------------------------------------------------------------
1527
1528 RooAbsPdf* pdf = (RooAbsPdf*) term->first() ;
1529
1530 RooAbsReal* partInt = pdf->createIntegral(termISet,termNSet,isetRangeName) ;
1531 partInt->setOperMode(operMode()) ;
1532 partInt->setStringAttribute("PROD_TERM_TYPE","IIIa") ;
1533
1534 isOwned=true ;
1535
1536 //cout << "processProductTerm(" << GetName() << ") case IIIa func = " << partInt->GetName() << endl ;
1537
1538 ret[0] = partInt ;
1539
1540 // Split mode results
1541 ret[1] = pdf->createIntegral(termISet,isetRangeName) ;
1542 ret[2] = pdf->createIntegral(termNSet,normRange()) ;
1543
1544 return ret ;
1545
1546
1547 } else {
1548
1549 // CASE IIIb: Normalized and partially integrated composite PDF term
1550 //---------------------------------------------------------------
1551
1552 // Use auxiliary class RooGenProdProj to calculate this term
1553 const std::string name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
1554 RooAbsReal* partInt = new RooGenProdProj(name.c_str(),name.c_str(),*term,termISet,termNSet,isetRangeName) ;
1555 partInt->setStringAttribute("PROD_TERM_TYPE","IIIb") ;
1556 partInt->setOperMode(operMode()) ;
1557
1558 //cout << "processProductTerm(" << GetName() << ") case IIIb func = " << partInt->GetName() << endl ;
1559
1560 isOwned=true ;
1561 ret[0] = partInt ;
1562
1563 const std::string name1 = makeRGPPName("PROD",*term,RooArgSet(),RooArgSet(),0) ;
1564
1565 // WVE FIX THIS
1566 RooProduct* tmp_prod = new RooProduct(name1.c_str(),name1.c_str(),*term) ;
1567
1568 ret[1] = tmp_prod->createIntegral(termISet,isetRangeName) ;
1569 ret[2] = tmp_prod->createIntegral(termNSet,normRange()) ;
1570
1571 return ret ;
1572 }
1573 }
1574
1575 // CASE IVa: Normalized non-integrated composite PDF term
1576 // -------------------------------------------------------
1577 if (nset && nset->getSize()>0 && term->getSize()>1) {
1578 // Composite term needs normalized integration
1579
1580 const std::string name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
1581 RooAbsReal* partInt = new RooGenProdProj(name.c_str(),name.c_str(),*term,termISet,termNSet,isetRangeName,normRange()) ;
1583
1584 partInt->setStringAttribute("PROD_TERM_TYPE","IVa") ;
1585 partInt->setOperMode(operMode()) ;
1586
1587 //cout << "processProductTerm(" << GetName() << ") case IVa func = " << partInt->GetName() << endl ;
1588
1589 isOwned=true ;
1590 ret[0] = partInt ;
1591
1592 const std::string name1 = makeRGPPName("PROD",*term,RooArgSet(),RooArgSet(),0) ;
1593
1594 // WVE FIX THIS
1595 RooProduct* tmp_prod = new RooProduct(name1.c_str(),name1.c_str(),*term) ;
1596
1597 ret[1] = tmp_prod->createIntegral(termISet,isetRangeName) ;
1598 ret[2] = tmp_prod->createIntegral(termNSet,normRange()) ;
1599
1600 return ret ;
1601 }
1602
1603 // CASE IVb: Normalized, non-integrated single PDF term
1604 // -----------------------------------------------------
1605 RooFIter pIter = term->fwdIterator() ;
1606 RooAbsPdf* pdf ;
1607 while((pdf=(RooAbsPdf*)pIter.next())) {
1608
1609 if (forceWrap) {
1610
1611 // Construct representative name of normalization wrapper
1612 TString name(pdf->GetName()) ;
1613 name.Append("_NORM[") ;
1614 RooFIter nIter = termNSet.fwdIterator() ;
1615 RooAbsArg* arg ;
1616 bool first(true) ;
1617 while((arg=(RooAbsArg*)nIter.next())) {
1618 if (!first) {
1619 name.Append(",") ;
1620 } else {
1621 first=false ;
1622 }
1623 name.Append(arg->GetName()) ;
1624 }
1625 if (normRange()) {
1626 name.Append("|") ;
1627 name.Append(normRange()) ;
1628 }
1629 name.Append("]") ;
1630
1631 RooAbsReal* partInt = new RooRealIntegral(name.Data(),name.Data(),*pdf,RooArgSet(),&termNSet) ;
1632 partInt->setStringAttribute("PROD_TERM_TYPE","IVb") ;
1633 isOwned=true ;
1634
1635 //cout << "processProductTerm(" << GetName() << ") case IVb func = " << partInt->GetName() << endl ;
1636
1637 ret[0] = partInt ;
1638
1639 ret[1] = pdf->createIntegral(RooArgSet()) ;
1640 ret[2] = pdf->createIntegral(termNSet,normRange()) ;
1641
1642 return ret ;
1643
1644
1645 } else {
1646 isOwned=false ;
1647
1648 //cout << "processProductTerm(" << GetName() << ") case IVb func = " << pdf->GetName() << endl ;
1649
1650
1651 pdf->setStringAttribute("PROD_TERM_TYPE","IVb") ;
1652 ret[0] = pdf ;
1653
1654 ret[1] = pdf->createIntegral(RooArgSet()) ;
1655 ret[2] = termNSet.getSize()>0 ? pdf->createIntegral(termNSet,normRange()) : ((RooAbsReal*)RooFit::RooConst(1).clone("1")) ;
1656 return ret ;
1657 }
1658 }
1659
1660 coutE(Eval) << "RooProdPdf::processProductTerm(" << GetName() << ") unidentified term!!!" << endl ;
1661 return ret ;
1662}
1663
1664
1665
1666
1667////////////////////////////////////////////////////////////////////////////////
1668/// Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
1669
1670std::string RooProdPdf::makeRGPPName(const char* pfx, const RooArgSet& term, const RooArgSet& iset,
1671 const RooArgSet& nset, const char* isetRangeName) const
1672{
1673 // Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
1674
1675 std::ostringstream os(pfx);
1676 os << "[";
1677
1678 RooFIter pIter = term.fwdIterator() ;
1679 // Encode component names
1680 bool first(true) ;
1681 RooAbsPdf* pdf ;
1682 while ((pdf=(RooAbsPdf*)pIter.next())) {
1683 if (!first) os << "_X_";
1684 first = false;
1685 os << pdf->GetName();
1686 }
1687 os << "]" << integralNameSuffix(iset,&nset,isetRangeName,true);
1688
1689 return os.str();
1690}
1691
1692
1693
1694////////////////////////////////////////////////////////////////////////////////
1695/// Force RooRealIntegral to offer all observables for internal integration
1696
1698{
1699 return true ;
1700}
1701
1702
1703
1704////////////////////////////////////////////////////////////////////////////////
1705/// Determine which part (if any) of given integral can be performed analytically.
1706/// If any analytical integration is possible, return integration scenario code.
1707///
1708/// RooProdPdf implements two strategies in implementing analytical integrals
1709///
1710/// First, PDF components whose entire set of dependents are requested to be integrated
1711/// can be dropped from the product, as they will integrate out to 1 by construction
1712///
1713/// Second, RooProdPdf queries each remaining component PDF for its analytical integration
1714/// capability of the requested set ('allVars'). It finds the largest common set of variables
1715/// that can be integrated by all remaining components. If such a set exists, it reconfirms that
1716/// each component is capable of analytically integrating the common set, and combines the components
1717/// individual integration codes into a single integration code valid for RooProdPdf.
1718
1720 const RooArgSet* normSet, const char* rangeName) const
1721{
1722 if (_forceNumInt) return 0 ;
1723
1724 // Declare that we can analytically integrate all requested observables
1725 analVars.add(allVars) ;
1726
1727 // Retrieve (or create) the required partial integral list
1728 Int_t code = getPartIntList(normSet,&allVars,rangeName);
1729
1730 return code+1 ;
1731}
1732
1733
1734
1735
1736////////////////////////////////////////////////////////////////////////////////
1737/// Return analytical integral defined by given scenario code
1738
1739double RooProdPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
1740{
1741 // No integration scenario
1742 if (code==0) {
1743 return getVal(normSet) ;
1744 }
1745
1746
1747 // WVE needs adaptation for rangename feature
1748
1749 // Partial integration scenarios
1750 CacheElem* cache = (CacheElem*) _cacheMgr.getObjByIndex(code-1) ;
1751
1752 // If cache has been sterilized, revive this slot
1753 if (cache==0) {
1754 std::unique_ptr<RooArgSet> vars{getParameters(RooArgSet())} ;
1755 RooArgSet nset = _cacheMgr.selectFromSet1(*vars, code-1) ;
1756 RooArgSet iset = _cacheMgr.selectFromSet2(*vars, code-1) ;
1757
1758 Int_t code2 = getPartIntList(&nset, &iset, rangeName) ;
1759
1760 // preceding call to getPartIntList guarantees non-null return
1761 // coverity[NULL_RETURNS]
1762 cache = (CacheElem*) _cacheMgr.getObj(&nset,&iset,&code2,rangeName) ;
1763 }
1764
1765 double val = calculate(*cache,true) ;
1766// cout << "RPP::aIWN(" << GetName() << ") ,code = " << code << ", value = " << val << endl ;
1767
1768 return val ;
1769}
1770
1771
1772
1773////////////////////////////////////////////////////////////////////////////////
1774/// Obsolete
1775
1776bool RooProdPdf::checkObservables(const RooArgSet* /*nset*/) const
1777{ return false ; }
1778
1779
1780
1781
1782////////////////////////////////////////////////////////////////////////////////
1783/// If this product contains exactly one extendable p.d.f return the extension abilities of
1784/// that p.d.f, otherwise return CanNotBeExtended
1785
1787{
1789}
1790
1791
1792
1793////////////////////////////////////////////////////////////////////////////////
1794/// Return the expected number of events associated with the extendable input PDF
1795/// in the product. If there is no extendable term, abort.
1796
1797double RooProdPdf::expectedEvents(const RooArgSet* nset) const
1798{
1799 if (_extendedIndex<0) {
1800 coutF(Generation) << "Requesting expected number of events from a RooProdPdf that does not contain an extended p.d.f" << endl ;
1801 throw std::logic_error(std::string("RooProdPdf ") + GetName() + " could not be extended.");
1802 }
1803
1804 return ((RooAbsPdf*)_pdfList.at(_extendedIndex))->expectedEvents(nset) ;
1805}
1806
1807
1808
1809////////////////////////////////////////////////////////////////////////////////
1810/// Return generator context optimized for generating events from product p.d.f.s
1811
1813 const RooArgSet* auxProto, bool verbose) const
1814{
1815 if (_useDefaultGen) return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ;
1816 return new RooProdGenContext(*this,vars,prototype,auxProto,verbose) ;
1817}
1818
1819
1820
1821////////////////////////////////////////////////////////////////////////////////
1822/// Query internal generation capabilities of component p.d.f.s and aggregate capabilities
1823/// into master configuration passed to the generator context
1824
1825Int_t RooProdPdf::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool staticInitOK) const
1826{
1827 if (!_useDefaultGen) return 0 ;
1828
1829 // Find the subset directVars that only depend on a single PDF in the product
1830 RooArgSet directSafe ;
1831 RooFIter dIter = directVars.fwdIterator() ;
1832 RooAbsArg* arg ;
1833 while((arg=(RooAbsArg*)dIter.next())) {
1834 if (isDirectGenSafe(*arg)) directSafe.add(*arg) ;
1835 }
1836
1837
1838 // Now find direct integrator for relevant components ;
1839 RooAbsPdf* pdf ;
1840 std::vector<Int_t> code;
1841 code.reserve(64);
1842 RooFIter pdfIter = _pdfList.fwdIterator();
1843 while((pdf=(RooAbsPdf*)pdfIter.next())) {
1844 RooArgSet pdfDirect ;
1845 Int_t pdfCode = pdf->getGenerator(directSafe,pdfDirect,staticInitOK);
1846 code.push_back(pdfCode);
1847 if (pdfCode != 0) {
1848 generateVars.add(pdfDirect) ;
1849 }
1850 }
1851
1852
1853 if (generateVars.getSize()>0) {
1854 Int_t masterCode = _genCode.store(code) ;
1855 return masterCode+1 ;
1856 } else {
1857 return 0 ;
1858 }
1859}
1860
1861
1862
1863////////////////////////////////////////////////////////////////////////////////
1864/// Forward one-time initialization call to component generation initialization
1865/// methods.
1866
1868{
1869 if (!_useDefaultGen) return ;
1870
1871 const std::vector<Int_t>& codeList = _genCode.retrieve(code-1) ;
1872 RooAbsPdf* pdf ;
1873 Int_t i(0) ;
1874 RooFIter pdfIter = _pdfList.fwdIterator();
1875 while((pdf=(RooAbsPdf*)pdfIter.next())) {
1876 if (codeList[i]!=0) {
1877 pdf->initGenerator(codeList[i]) ;
1878 }
1879 i++ ;
1880 }
1881}
1882
1883
1884
1885////////////////////////////////////////////////////////////////////////////////
1886/// Generate a single event with configuration specified by 'code'
1887/// Defer internal generation to components as encoded in the _genCode
1888/// registry for given generator code.
1889
1891{
1892 if (!_useDefaultGen) return ;
1893
1894 const std::vector<Int_t>& codeList = _genCode.retrieve(code-1) ;
1895 RooAbsPdf* pdf ;
1896 Int_t i(0) ;
1897 RooFIter pdfIter = _pdfList.fwdIterator();
1898 while((pdf=(RooAbsPdf*)pdfIter.next())) {
1899 if (codeList[i]!=0) {
1900 pdf->generateEvent(codeList[i]) ;
1901 }
1902 i++ ;
1903 }
1904}
1905
1906
1907
1908////////////////////////////////////////////////////////////////////////////////
1909/// Return RooAbsArg components contained in the cache
1910
1912{
1913 RooArgList ret ;
1914 ret.add(_partList) ;
1915 ret.add(_numList) ;
1916 ret.add(_denList) ;
1917 if (_rearrangedNum) ret.add(*_rearrangedNum) ;
1918 if (_rearrangedDen) ret.add(*_rearrangedDen) ;
1919 return ret ;
1920
1921}
1922
1923
1924
1925////////////////////////////////////////////////////////////////////////////////
1926/// Hook function to print cache contents in tree printing of RooProdPdf
1927
1928void RooProdPdf::CacheElem::printCompactTreeHook(ostream& os, const char* indent, Int_t curElem, Int_t maxElem)
1929{
1930 if (curElem==0) {
1931 os << indent << "RooProdPdf begin partial integral cache" << endl ;
1932 }
1933
1934 auto indent2 = std::string(indent) + "[" + std::to_string(curElem) + "]";
1935 for(auto const& arg : _partList) {
1936 arg->printCompactTree(os,indent2.c_str()) ;
1937 }
1938
1939 if (curElem==maxElem) {
1940 os << indent << "RooProdPdf end partial integral cache" << endl ;
1941 }
1942}
1943
1944
1945
1946////////////////////////////////////////////////////////////////////////////////
1947/// Forward determination of safety of internal generator code to
1948/// component p.d.f that would generate the given observable
1949
1951{
1952 // Only override base class behaviour if default generator method is enabled
1953 if (!_useDefaultGen) return RooAbsPdf::isDirectGenSafe(arg) ;
1954
1955 // Argument may appear in only one PDF component
1956 RooAbsPdf* pdf, *thePdf(0) ;
1957 RooFIter pdfIter = _pdfList.fwdIterator();
1958 while((pdf=(RooAbsPdf*)pdfIter.next())) {
1959
1960 if (pdf->dependsOn(arg)) {
1961 // Found PDF depending on arg
1962
1963 // If multiple PDFs depend on arg directGen is not safe
1964 if (thePdf) return false ;
1965
1966 thePdf = pdf ;
1967 }
1968 }
1969 // Forward call to relevant component PDF
1970 return thePdf?(thePdf->isDirectGenSafe(arg)):false ;
1971}
1972
1973
1974
1975////////////////////////////////////////////////////////////////////////////////
1976/// Look up user specified normalization set for given input PDF component
1977
1979{
1980 Int_t idx = _pdfList.index(&pdf) ;
1981 if (idx<0) return nullptr;
1982 return _pdfNSetList[idx].get() ;
1983}
1984
1985
1986
1987////////////////////////////////////////////////////////////////////////////////
1988/// Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
1989/// The observables set is required to distinguish unambiguously p.d.f in terms
1990/// of observables and parameters, which are not constraints, and p.d.fs in terms
1991/// of parameters only, which can serve as constraints p.d.f.s
1992
1993RooArgSet* RooProdPdf::getConstraints(const RooArgSet& observables, RooArgSet& constrainedParams, bool stripDisconnected) const
1994{
1995 RooArgSet constraints ;
1996 RooArgSet pdfParams, conParams ;
1997
1998 // Loop over p.d.f. components
1999 for(auto * pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
2000 // A constraint term is a p.d.f that does not depend on any of the listed observables
2001 // but does depends on any of the parameters that should be constrained
2002 RooArgSet tmp;
2003 if (!pdf->dependsOnValue(observables) && pdf->dependsOnValue(constrainedParams)) {
2004 constraints.add(*pdf) ;
2005 pdf->getParameters(&observables, tmp);
2006 conParams.add(tmp,true) ;
2007 } else {
2008 pdf->getParameters(&observables, tmp);
2009 pdfParams.add(tmp,true) ;
2010 }
2011 }
2012
2013 // Strip any constraints that are completely decoupled from the other product terms
2014 RooArgSet* finalConstraints = new RooArgSet("constraints") ;
2015 for(auto * pdf : static_range_cast<RooAbsPdf*>(constraints)) {
2016 if (pdf->dependsOnValue(pdfParams) || !stripDisconnected) {
2017 finalConstraints->add(*pdf) ;
2018 } else {
2019 coutI(Minimization) << "RooProdPdf::getConstraints(" << GetName() << ") omitting term " << pdf->GetName()
2020 << " as constraint term as it does not share any parameters with the other pdfs in product. "
2021 << "To force inclusion in likelihood, add an explicit Constrain() argument for the target parameter" << endl ;
2022 }
2023 }
2024
2025 // Now remove from constrainedParams all parameters that occur exclusively in constraint term and not in regular pdf term
2026
2027 std::unique_ptr<RooArgSet> cexl{static_cast<RooArgSet*>(conParams.selectCommon(constrainedParams))};
2028 cexl->remove(pdfParams,true,true) ;
2029 constrainedParams.remove(*cexl,true,true) ;
2030
2031 return finalConstraints ;
2032}
2033
2034
2035
2036
2037////////////////////////////////////////////////////////////////////////////////
2038/// Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
2039/// The observables set is required to distinguish unambiguously p.d.f in terms
2040/// of observables and parameters, which are not constraints, and p.d.fs in terms
2041/// of parameters only, which can serve as constraints p.d.f.s
2042
2044{
2045 RooArgSet* connectedPars = new RooArgSet("connectedPars") ;
2046 for (const auto arg : _pdfList) {
2047 // Check if term is relevant
2048 if (arg->dependsOn(observables)) {
2049 RooArgSet tmp;
2050 arg->getParameters(&observables, tmp);
2051 connectedPars->add(tmp) ;
2052 }
2053 }
2054 return connectedPars ;
2055}
2056
2057
2058
2059
2060////////////////////////////////////////////////////////////////////////////////
2061
2062void RooProdPdf::getParametersHook(const RooArgSet* nset, RooArgSet* params, bool stripDisconnected) const
2063{
2064 if (!stripDisconnected) return ;
2065 if (!nset || nset->getSize()==0) return ;
2066
2067 // Get/create appropriate term list for this normalization set
2068 Int_t code = getPartIntList(nset, nullptr);
2069 RooArgList & plist = static_cast<CacheElem*>(_cacheMgr.getObj(nset, &code))->_partList;
2070
2071 // Strip any terms from params that do not depend on any term
2072 RooArgSet tostrip ;
2073 for (auto param : *params) {
2074 bool anyDep(false) ;
2075 for (auto term : plist) {
2076 if (term->dependsOnValue(*param)) {
2077 anyDep=true ;
2078 }
2079 }
2080 if (!anyDep) {
2081 tostrip.add(*param) ;
2082 }
2083 }
2084
2085 if (tostrip.getSize()>0) {
2086 params->remove(tostrip,true,true);
2087 }
2088
2089}
2090
2091
2092
2093////////////////////////////////////////////////////////////////////////////////
2094/// Interface function used by test statistics to freeze choice of range
2095/// for interpretation of conditional product terms
2096
2097void RooProdPdf::selectNormalizationRange(const char* rangeName, bool force)
2098{
2099 if (!force && _refRangeName) {
2100 return ;
2101 }
2102
2103 fixRefRange(rangeName) ;
2104}
2105
2106
2107
2108
2109////////////////////////////////////////////////////////////////////////////////
2110
2111void RooProdPdf::fixRefRange(const char* rangeName)
2112{
2113 _refRangeName = (TNamed*)RooNameReg::ptr(rangeName) ;
2114}
2115
2116
2117
2118////////////////////////////////////////////////////////////////////////////////
2119/// Forward the plot sampling hint from the p.d.f. that defines the observable obs
2120
2121std::list<double>* RooProdPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
2122{
2123 RooAbsPdf* pdf ;
2124 RooFIter pdfIter = _pdfList.fwdIterator();
2125 while((pdf=(RooAbsPdf*)pdfIter.next())) {
2126 list<double>* hint = pdf->plotSamplingHint(obs,xlo,xhi) ;
2127 if (hint) {
2128 return hint ;
2129 }
2130 }
2131
2132 return nullptr;
2133}
2134
2135
2136
2137////////////////////////////////////////////////////////////////////////////////
2138/// If all components that depend on obs are binned that so is the product
2139
2141{
2142 RooAbsPdf* pdf ;
2143 RooFIter pdfIter = _pdfList.fwdIterator();
2144 while((pdf=(RooAbsPdf*)pdfIter.next())) {
2145 if (pdf->dependsOn(obs) && !pdf->isBinnedDistribution(obs)) {
2146 return false ;
2147 }
2148 }
2149
2150 return true ;
2151}
2152
2153
2154
2155
2156
2157
2158////////////////////////////////////////////////////////////////////////////////
2159/// Forward the plot sampling hint from the p.d.f. that defines the observable obs
2160
2161std::list<double>* RooProdPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
2162{
2163 RooAbsPdf* pdf ;
2164 RooFIter pdfIter = _pdfList.fwdIterator();
2165 while((pdf=(RooAbsPdf*)pdfIter.next())) {
2166 list<double>* hint = pdf->binBoundaries(obs,xlo,xhi) ;
2167 if (hint) {
2168 return hint ;
2169 }
2170 }
2171
2172 return nullptr;
2173}
2174
2175
2176////////////////////////////////////////////////////////////////////////////////
2177/// Label OK'ed components of a RooProdPdf with cache-and-track, _and_ label all RooProdPdf
2178/// descendants with extra information about (conditional) normalization, needed to be able
2179/// to Cache-And-Track them outside the RooProdPdf context.
2180
2182{
2183 for (const auto parg : _pdfList) {
2184
2185 if (parg->canNodeBeCached()==Always) {
2186 trackNodes.add(*parg) ;
2187// cout << "tracking node RooProdPdf component " << parg << " " << parg->IsA()->GetName() << "::" << parg->GetName() << endl ;
2188
2189 // Additional processing to fix normalization sets in case product defines conditional observables
2190 RooArgSet* pdf_nset = findPdfNSet((RooAbsPdf&)(*parg)) ;
2191 if (pdf_nset) {
2192 // Check if conditional normalization is specified
2194 if (string("nset")==pdf_nset->GetName() && pdf_nset->getSize()>0) {
2195 parg->setStringAttribute("CATNormSet",getColonSeparatedNameString(*pdf_nset).c_str()) ;
2196 }
2197 if (string("cset")==pdf_nset->GetName()) {
2198 parg->setStringAttribute("CATCondSet",getColonSeparatedNameString(*pdf_nset).c_str()) ;
2199 }
2200 } else {
2201 coutW(Optimization) << "RooProdPdf::setCacheAndTrackHints(" << GetName() << ") WARNING product pdf does not specify a normalization set for component " << parg->GetName() << endl ;
2202 }
2203 }
2204 }
2205}
2206
2207
2208
2209////////////////////////////////////////////////////////////////////////////////
2210/// Customized printing of arguments of a RooProdPdf to more intuitively reflect the contents of the
2211/// product operator construction
2212
2213void RooProdPdf::printMetaArgs(ostream& os) const
2214{
2215 for (int i=0 ; i<_pdfList.getSize() ; i++) {
2216 if (i>0) os << " * " ;
2217 RooArgSet* ncset = _pdfNSetList[i].get() ;
2218 os << _pdfList.at(i)->GetName() ;
2219 if (ncset->getSize()>0) {
2220 if (string("nset")==ncset->GetName()) {
2221 os << *ncset ;
2222 } else {
2223 os << "|" ;
2224 RooFIter nciter = ncset->fwdIterator() ;
2225 RooAbsArg* arg ;
2226 bool first(true) ;
2227 while((arg=(RooAbsArg*)nciter.next())) {
2228 if (!first) {
2229 os << "," ;
2230 } else {
2231 first = false ;
2232 }
2233 os << arg->GetName() ;
2234 }
2235 }
2236 }
2237 }
2238 os << " " ;
2239}
2240
2241
2242
2243////////////////////////////////////////////////////////////////////////////////
2244/// Implement support for node removal
2245
2246bool RooProdPdf::redirectServersHook(const RooAbsCollection& newServerList, bool /*mustReplaceAll*/, bool nameChange, bool /*isRecursive*/)
2247{
2248 if (nameChange && _pdfList.find("REMOVAL_DUMMY")) {
2249
2250 cxcoutD(LinkStateMgmt) << "RooProdPdf::redirectServersHook(" << GetName() << "): removing REMOVAL_DUMMY" << endl ;
2251
2252 // Remove node from _pdfList proxy and remove corresponding entry from normset list
2253 RooAbsArg* pdfDel = _pdfList.find("REMOVAL_DUMMY") ;
2254
2255 _pdfNSetList.erase(_pdfNSetList.begin() + _pdfList.index("REMOVAL_DUMMY")) ;
2256 _pdfList.remove(*pdfDel) ;
2257
2258 // Clear caches
2259 _cacheMgr.reset() ;
2260 }
2261
2262 // If the replaced server is an observable that is used in any of the
2263 // normalization sets for conditional fits, replace the element in the
2264 // normalization set too.
2265 for(std::unique_ptr<RooArgSet> const& normSet : _pdfNSetList) {
2266 for(RooAbsArg * arg : *normSet) {
2267 if(RooAbsArg * newArg = arg->findNewServer(newServerList, nameChange)) {
2268 // Need to do some tricks here because it's not possible to replace in
2269 // an owning RooAbsCollection.
2270 normSet->releaseOwnership();
2271 normSet->replace(*arg, *newArg->cloneTree());
2272 normSet->takeOwnership();
2273 delete arg;
2274 }
2275 }
2276 }
2277
2278 return false ;
2279}
2280
2281void RooProdPdf::CacheElem::writeToStream(std::ostream& os) const {
2282 using namespace RooHelpers;
2283 os << "_partList\n";
2284 os << getColonSeparatedNameString(_partList) << "\n";
2285 os << "_numList\n";
2286 os << getColonSeparatedNameString(_numList) << "\n";
2287 os << "_denList\n";
2288 os << getColonSeparatedNameString(_denList) << "\n";
2289 os << "_ownedList\n";
2290 os << getColonSeparatedNameString(_ownedList) << "\n";
2291 os << "_normList\n";
2292 for(auto const& set : _normList) {
2293 os << getColonSeparatedNameString(*set) << "\n";
2294 }
2295 os << "_isRearranged" << "\n";
2296 os << _isRearranged << "\n";
2297 os << "_rearrangedNum" << "\n";
2298 if(_rearrangedNum) {
2299 os << getColonSeparatedNameString(*_rearrangedNum) << "\n";
2300 } else {
2301 os << "nullptr" << "\n";
2302 }
2303 os << "_rearrangedDen" << "\n";
2304 if(_rearrangedDen) {
2305 os << getColonSeparatedNameString(*_rearrangedDen) << "\n";
2306 } else {
2307 os << "nullptr" << "\n";
2308 }
2309}
2310
2311void RooProdPdf::writeCacheToStream(std::ostream& os, RooArgSet const* nset) const {
2312 getCacheElem(nset)->writeToStream(os);
2313}
2314
2315void RooProdPdf::fillNormSetForServer(RooArgSet const& normSet, RooAbsArg const& server, RooArgSet& serverNormSet) const {
2316 auto * pdfNset = findPdfNSet(static_cast<RooAbsPdf const&>(server));
2317 if (pdfNset && !pdfNset->empty()) {
2318 for(auto * arg : *pdfNset) if(server.dependsOn(*arg)) serverNormSet.add(*arg);
2319 } else {
2320 for(auto * arg : normSet) if(server.dependsOn(*arg)) serverNormSet.add(*arg);
2321 }
2322}
#define coutI(a)
Definition: RooMsgService.h:34
#define cxcoutD(a)
Definition: RooMsgService.h:85
#define coutW(a)
Definition: RooMsgService.h:36
#define dologD(a)
Definition: RooMsgService.h:69
#define coutF(a)
Definition: RooMsgService.h:38
#define coutE(a)
Definition: RooMsgService.h:37
#define TRACE_DESTROY
Definition: RooTrace.h:24
#define TRACE_CREATE
Definition: RooTrace.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
static void indent(ostringstream &buf, int indent_level)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
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
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2447
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...
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:77
RooExpensiveObjectCache & expensiveObjectCache() const
Definition: RooAbsArg.cxx:2354
TObject * Clone(const char *newname=0) const override
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:89
void setOperMode(OperMode mode, bool recurseADirty=true)
Set the operation mode of this node.
Definition: RooAbsArg.cxx:1876
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
Definition: RooAbsArg.cxx:314
void treeNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=0, bool doBranch=true, bool doLeaf=true, bool valueOnly=false, bool recurseNonDerived=false) const
Fill supplied list with nodes of the arg tree, following all server links, starting with ourself as t...
Definition: RooAbsArg.cxx:527
friend class RooArgSet
Definition: RooAbsArg.h:648
virtual void setExpensiveObjectCache(RooExpensiveObjectCache &cache)
Definition: RooAbsArg.h:526
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2275
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
Definition: RooAbsArg.cxx:327
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:835
void setAttribute(const Text_t *name, bool value=true)
Set (default) or clear a named boolean attribute of this object.
Definition: RooAbsArg.cxx:282
RooArgSet * getComponents() const
Create a RooArgSet with all components (branch nodes) of the expression tree headed by this object.
Definition: RooAbsArg.cxx:782
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:569
RooAbsArg * findNewServer(const RooAbsCollection &newSet, bool nameChange) const
Find the new server in the specified set that matches the old server.
Definition: RooAbsArg.cxx:1150
OperMode operMode() const
Query the operation mode of this node.
Definition: RooAbsArg.h:507
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
Storage_t const & get() const
Const access to the underlying stl container.
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
void printValue(std::ostream &os) const override
Print value of collection, i.e.
Int_t getSize() const
Return the number of elements in the collection.
const char * GetName() const override
Returns name of object.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
RooFIter fwdIterator() const
One-time forward iterator.
bool overlaps(Iterator_t otherCollBegin, Iterator_t otherCollEnd) const
Storage_t::size_type size() const
RooAbsArg * first() const
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
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.
friend class CacheElem
The cache manager.
Definition: RooAbsPdf.h:366
virtual void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
Definition: RooAbsPdf.cxx:2324
TString _normRange
Normalization range.
Definition: RooAbsPdf.h:392
virtual bool isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
Definition: RooAbsPdf.cxx:2337
RooArgSet const * _normSet
Normalization integral (owned by _normMgr)
Definition: RooAbsPdf.h:354
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:263
friend class RooRealIntegral
Definition: RooAbsPdf.h:346
@ CanNotBeExtended
Definition: RooAbsPdf.h:257
const char * normRange() const
Definition: RooAbsPdf.h:292
virtual void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
Definition: RooAbsPdf.cxx:2312
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, bool verbose=false) const
Interface function to create a generator context from a p.d.f.
Definition: RooAbsPdf.cxx:1903
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooAbsPdf.cxx:2302
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
virtual double getMax(const char *name=0) const
Get maximum of currently defined range.
virtual double getMin(const char *name=0) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
virtual std::list< double > * binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
Retrieve bin boundaries if this distribution is binned in obs.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:94
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:553
bool _forceNumInt
Force numerical integration if flag set.
Definition: RooAbsReal.h:487
TString integralNameSuffix(const RooArgSet &iset, const RooArgSet *nset=0, const char *rangeName=0, bool omitEmpty=false) const
Construct string with unique suffix name to give to integral object that encodes integrated observabl...
Definition: RooAbsReal.cxx:797
virtual bool isBinnedDistribution(const RooArgSet &) const
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
Definition: RooAbsReal.h:345
virtual std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
Definition: RooAddition.h:27
const RooArgList & list1() const
Definition: RooAddition.h:42
static TClass * Class()
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:57
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:180
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Getter function without integration set.
RooArgSet selectFromSet1(RooArgSet const &argSet, int index) const
Create RooArgSet contatining the objects that are both in the cached set 1.
T * getObjByIndex(Int_t index) const
Retrieve payload object by slot index.
RooArgSet selectFromSet2(RooArgSet const &argSet, int index) const
Create RooArgSet contatining the objects that are both in the cached set 2.
void reset()
Clear the cache.
Int_t lastIndex() const
Return index of slot used in last get or set operation.
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
Setter function without integration set.
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:26
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...
bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false) override
Remove object 'var' from set and deregister 'var' as server to owner.
TObject * clone(const char *newname) const override
Definition: RooConstVar.h:32
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
Definition: RooCustomizer.h:35
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition: RooFormulaVar.h:30
RooAbsArg * getParameter(const char *name) const
Return pointer to parameter with given name.
Definition: RooFormulaVar.h:44
RooGenProdProj is an auxiliary class for RooProdPdf that calculates a general normalised projection o...
RooLinkedList is an collection class for internal use, storing a collection of RooAbsArg pointers in ...
Definition: RooLinkedList.h:38
TObject * At(int index) const
Return object stored in sequential position given by index.
RooLinkedListIterImpl end() const
void Delete(Option_t *o=0) override
Remove all elements in collection and delete all elements NB: Collection does not own elements,...
virtual void Add(TObject *arg)
Definition: RooLinkedList.h:67
RooLinkedListIterImpl begin() const
Int_t IndexOf(const char *name) const
Return position of given object in list.
static const char * str(const TNamed *ptr)
Return C++ string corresponding to given TNamed pointer.
Definition: RooNameReg.cxx:102
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:92
RooArgList containedArgs(Action) override
Return RooAbsArg components contained in the cache.
std::unique_ptr< RooAbsReal > _rearrangedNum
Definition: RooProdPdf.h:147
void printCompactTreeHook(std::ostream &, const char *, Int_t, Int_t) override
Hook function to print cache contents in tree printing of RooProdPdf.
RooArgList _ownedList
Definition: RooProdPdf.h:144
RooArgList _partList
Definition: RooProdPdf.h:141
std::vector< std::unique_ptr< RooArgSet > > _normList
Definition: RooProdPdf.h:145
std::unique_ptr< RooAbsReal > _rearrangedDen
Definition: RooProdPdf.h:148
void writeToStream(std::ostream &os) const
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:33
void setCacheAndTrackHints(RooArgSet &) override
Label OK'ed components of a RooProdPdf with cache-and-track, and label all RooProdPdf descendants wit...
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Query internal generation capabilities of component p.d.f.s and aggregate capabilities into master co...
void rearrangeProduct(CacheElem &) const
void factorizeProduct(const RooArgSet &normSet, const RooArgSet &intSet, RooLinkedList &termList, RooLinkedList &normList, RooLinkedList &impDepList, RooLinkedList &crossDepList, RooLinkedList &intList) const
Factorize product in irreducible terms for given choice of integration/normalization.
Definition: RooProdPdf.cxx:548
~RooProdPdf() override
Destructor.
Definition: RooProdPdf.cxx:427
Int_t _extendedIndex
Index of extended PDF (if any)
Definition: RooProdPdf.h:173
Int_t getPartIntList(const RooArgSet *nset, const RooArgSet *iset, const char *isetRangeName=0) const
Return list of (partial) integrals of product terms for integration of p.d.f over observables iset wh...
Definition: RooProdPdf.cxx:732
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Evaluate product of PDFs in batch mode.
Definition: RooProdPdf.cxx:491
RooProdPdf()
Default constructor.
Definition: RooProdPdf.cxx:82
RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, bool verbose=false) const override
Return generator context optimized for generating events from product p.d.f.s.
bool checkObservables(const RooArgSet *nset) const override
Obsolete.
bool _useDefaultGen
Use default or distributed event generator.
Definition: RooProdPdf.h:176
std::vector< std::unique_ptr< RooArgSet > > _pdfNSetList
List of PDF component normalization sets.
Definition: RooProdPdf.h:172
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=0) const override
Determine which part (if any) of given integral can be performed analytically.
bool forceAnalyticalInt(const RooAbsArg &dep) const override
Force RooRealIntegral to offer all observables for internal integration.
void writeCacheToStream(std::ostream &os, RooArgSet const *nset) const
void getParametersHook(const RooArgSet *, RooArgSet *, bool stripDisconnected) const override
TNamed * _refRangeName
Reference range name for interpretation of conditional products.
Definition: RooProdPdf.h:178
RooAICRegistry _genCode
! Registry of composite direct generator codes
Definition: RooProdPdf.h:168
RooListProxy _pdfList
List of PDF components.
Definition: RooProdPdf.h:171
void selectNormalizationRange(const char *rangeName=0, bool force=false) override
Interface function used by test statistics to freeze choice of range for interpretation of conditiona...
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooProdPdf to more intuitively reflect the contents of the prod...
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Forward the plot sampling hint from the p.d.f. that defines the observable obs.
RooAbsReal * specializeIntegral(RooAbsReal &orig, const char *targetRangeName) const
void fillNormSetForServer(RooArgSet const &normSet, RooAbsArg const &server, RooArgSet &serverNormSet) const override
Fills a RooArgSet to be used as the normalization set for a server, given a normalization set for thi...
RooObjCacheManager _cacheMgr
Definition: RooProdPdf.h:154
std::string makeRGPPName(const char *pfx, const RooArgSet &term, const RooArgSet &iset, const RooArgSet &nset, const char *isetRangeName) const
Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
bool isDirectGenSafe(const RooAbsArg &arg) const override
Forward determination of safety of internal generator code to component p.d.f that would generate the...
RooArgSet * findPdfNSet(RooAbsPdf const &pdf) const
Look up user specified normalization set for given input PDF component.
RooAbsReal * makeCondPdfRatioCorr(RooAbsReal &term, const RooArgSet &termNset, const RooArgSet &termImpSet, const char *normRange, const char *refRange) const
For single normalization ranges.
ExtendMode extendMode() const override
If this product contains exactly one extendable p.d.f return the extension abilities of that p....
std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override
Forward the plot sampling hint from the p.d.f. that defines the observable obs.
double expectedEvents(const RooArgSet *nset) const override
Return the expected number of events associated with the extendable input PDF in the product.
std::vector< RooAbsReal * > processProductTerm(const RooArgSet *nset, const RooArgSet *iset, const char *isetRangeName, const RooArgSet *term, const RooArgSet &termNSet, const RooArgSet &termISet, bool &isOwned, bool forceWrap=false) const
Calculate integrals of factorized product terms over observables iset while normalized to observables...
TObject * clone(const char *newname) const override
Definition: RooProdPdf.h:56
bool isBinnedDistribution(const RooArgSet &obs) const override
If all components that depend on obs are binned that so is the product.
friend class RooProdGenContext
Definition: RooProdPdf.h:163
RooArgSet * getConnectedParameters(const RooArgSet &observables) const
Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
double calculate(const RooProdPdf::CacheElem &cache, bool verbose=false) const
Calculate running product of pdfs terms, using the supplied normalization set in 'normSetList' for ea...
Definition: RooProdPdf.cxx:460
RooArgSet _defNormSet
Default normalization set.
Definition: RooProdPdf.h:181
RooArgSet * getConstraints(const RooArgSet &observables, RooArgSet &constrainedParams, bool stripDisconnected) const override
Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
RooAbsReal * specializeRatio(RooFormulaVar &input, const char *targetRangeName) const
CacheElem * getCacheElem(RooArgSet const *nset) const
The cache manager.
Definition: RooProdPdf.cxx:433
bool redirectServersHook(const RooAbsCollection &, bool, bool, bool) override
Implement support for node removal.
void initGenerator(Int_t code) override
Forward one-time initialization call to component generation initialization methods.
void generateEvent(Int_t code) override
Generate a single event with configuration specified by 'code' Defer internal generation to component...
void groupProductTerms(std::list< std::vector< RooArgSet * > > &groupedTerms, RooArgSet &outerIntDeps, const RooLinkedList &terms, const RooLinkedList &norms, const RooLinkedList &imps, const RooLinkedList &ints, const RooLinkedList &cross) const
Group product into terms that can be calculated independently.
void fixRefRange(const char *rangeName)
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const override
Return analytical integral defined by given scenario code.
double evaluate() const override
Calculate current value of object.
Definition: RooProdPdf.cxx:449
void initializeFromCmdArgList(const RooArgSet &fullPdfSet, const RooLinkedList &l)
Initialize RooProdPdf configuration from given list of RooCmdArg configuration arguments and set of '...
Definition: RooProdPdf.cxx:368
double _cutOff
Cutoff parameter for running product.
Definition: RooProdPdf.h:170
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:29
static TClass * Class()
RooArgList components()
Definition: RooProduct.h:47
RooRealIntegral performs hybrid numerical/analytical integrals of RooAbsReal objects.
RooArgSet intVars() const
const RooAbsReal & integrand() const
static TClass * Class()
const char * intRange()
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
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
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:445
Basic string class.
Definition: TString.h:136
Ssiz_t Length() const
Definition: TString.h:410
const char * Data() const
Definition: TString.h:369
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
RooConstVar & RooConst(double val)
Double_t x[n]
Definition: legend1.C:17
double T(double x)
Definition: ChebyshevPol.h:34
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
std::map< DataKey, RooSpan< const double > > DataMap
Definition: DataMap.h:59
@ Minimization
Definition: RooGlobalFunc.h:63
@ Generation
Definition: RooGlobalFunc.h:63
@ Optimization
Definition: RooGlobalFunc.h:64
@ InputArguments
Definition: RooGlobalFunc.h:64
@ LinkStateMgmt
Definition: RooGlobalFunc.h:63
std::string getColonSeparatedNameString(RooArgSet const &argSet)
Create a string with all sorted names of RooArgSet elements separated by colons.
Definition: RooHelpers.cxx:252
Definition: first.py:1
auto * l
Definition: textangle.C:4
static void output()