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