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