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