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
412////////////////////////////////////////////////////////////////////////////////
413/// Evaluate product of PDFs in batch mode.
415{
416 if (cache._isRearranged) {
417 auto numerator = ctx.at(cache._rearrangedNum.get());
418 auto denominator = ctx.at(cache._rearrangedDen.get());
419 RooBatchCompute::compute(ctx.config(caller), RooBatchCompute::Ratio, ctx.output(), {numerator, denominator});
420 } else {
421 std::vector<std::span<const double>> factors;
422 factors.reserve(cache._partList.size());
423 for (const RooAbsArg *i : cache._partList) {
424 auto span = ctx.at(i);
425 factors.push_back(span);
426 }
427 std::array<double, 1> special{static_cast<double>(factors.size())};
429 }
430}
431
432namespace {
433
434template<class T>
435void eraseNullptrs(std::vector<T*>& v) {
436 v.erase(std::remove_if(v.begin(), v.end(), [](T* x){ return x == nullptr; } ), v.end());
437}
438
439void removeCommon(std::vector<RooAbsArg*> &v, std::span<RooAbsArg * const> other) {
440
441 for (auto const& arg : other) {
442 auto namePtrMatch = [&arg](const RooAbsArg* elm) {
443 return elm != nullptr && elm->namePtr() == arg->namePtr();
444 };
445
446 auto found = std::find_if(v.begin(), v.end(), namePtrMatch);
447 if(found != v.end()) {
448 *found = nullptr;
449 }
450 }
452}
453
454void addCommon(std::vector<RooAbsArg*> &v, std::vector<RooAbsArg*> const& o1, std::vector<RooAbsArg*> const& o2) {
455
456 for (auto const& arg : o1) {
457 auto namePtrMatch = [&arg](const RooAbsArg* elm) {
458 return elm->namePtr() == arg->namePtr();
459 };
460
461 if(std::find_if(o2.begin(), o2.end(), namePtrMatch) != o2.end()) {
462 v.push_back(arg);
463 }
464 }
465}
466
467bool isRangeIdentical(RooArgSet const &observables, TString const &normRange, TNamed *refRangeName)
468{
469 // FK: Here the refRange should be compared to normRange, if it's set, and to the normObs range if it's not set
470 const char *range = normRange.Length() > 0 ? normRange.Data() : nullptr;
471 const char *refRange = RooNameReg::str(refRangeName);
472 for (auto const *normObs : static_range_cast<RooRealVar *>(observables)) {
473 if (normObs->getMin(range) != normObs->getMin(refRange) || normObs->getMax(range) != normObs->getMax(refRange))
474 return false;
475 }
476 return true;
477}
478
479}
480
481
482////////////////////////////////////////////////////////////////////////////////
483/// Factorize product in irreducible terms for given choice of integration/normalization
484
486{
487 // List of all term dependents: normalization and imported
488 std::vector<RooArgSet> depAllList;
489 std::vector<RooArgSet> depIntNoNormList;
490
491 // Setup lists for factorization terms and their dependents
492 RooArgSet* term(nullptr);
493 RooArgSet* termIntDeps(nullptr);
495
496 std::vector<RooAbsArg*> pdfIntNoNormDeps;
497 std::vector<RooAbsArg*> pdfIntSet;
498 std::vector<RooAbsArg*> pdfNSet;
499 std::vector<RooAbsArg*> pdfCSet;
500 std::vector<RooAbsArg*> pdfNormDeps; // Dependents to be normalized for the PDF
501 std::vector<RooAbsArg*> pdfAllDeps; // All dependents of this PDF
502
503 // Loop over the PDFs
504 for(std::size_t iPdf = 0; iPdf < _pdfList.size(); ++iPdf) {
505 RooAbsPdf& pdf = static_cast<RooAbsPdf&>(_pdfList[iPdf]);
507
508 pdfNSet.clear();
509 pdfCSet.clear();
510
511 // Make iterator over tree leaf node list to get the observables.
512 // This code is borrowed from RooAbsPdf::getObservables().
513 // RooAbsArg::treeNodeServer list is relatively expensive, so we only do it
514 // once and use it in a lambda function.
515 RooArgSet pdfLeafList("leafNodeServerList") ;
516 pdf.treeNodeServerList(&pdfLeafList,nullptr,false,true,true) ;
518 std::vector<RooAbsArg*> & out,
519 const RooArgSet& dataList) {
520 for (const auto arg : pdfLeafList) {
521 if (arg->dependsOnValue(dataList) && arg->isLValue()) {
522 out.push_back(arg) ;
523 }
524 }
525 };
526
527 // Reduce pdfNSet to actual dependents
528 if (0 == strcmp("cset", pdfNSetOrig.GetName())) {
531 pdfCSet = pdfNSetOrig.get();
532 } else {
533 // Interpret at NSet
535 }
536
537
538 pdfNormDeps.clear();
539 pdfAllDeps.clear();
540
541 // Make list of all dependents of this PDF
543
544
545// std::cout << GetName() << ": pdf = " << pdf->GetName() << " pdfAllDeps = " << pdfAllDeps << " pdfNSet = " << *pdfNSet << " pdfCSet = " << *pdfCSet << std::endl;
546
547 // Make list of normalization dependents for this PDF;
548 if (!pdfNSet.empty()) {
549 // PDF is conditional
551 } else {
552 // PDF is regular
554 }
555
556// std::cout << GetName() << ": pdfNormDeps for " << pdf->GetName() << " = " << pdfNormDeps << std::endl;
557
558 pdfIntSet.clear();
560
561 // WVE if we have no norm deps, conditional observables should be taken out of pdfIntSet
562 if (pdfNormDeps.empty() && !pdfCSet.empty()) {
564// std::cout << GetName() << ": have no norm deps, removing conditional observables from intset" << std::endl;
565 }
566
567 pdfIntNoNormDeps.clear();
570
571// std::cout << GetName() << ": pdf = " << pdf->GetName() << " intset = " << *pdfIntSet << " pdfIntNoNormDeps = " << pdfIntNoNormDeps << std::endl;
572
573 // Check if this PDF has dependents overlapping with one of the existing terms
574 bool done = false;
575 int j = 0;
576 auto lIter = factorized.terms.begin();
577 auto ldIter = factorized.norms.begin();
578 for(;lIter != factorized.terms.end(); (++lIter, ++ldIter, ++j)) {
579 RooArgSet *termNormDeps = static_cast<RooArgSet*>(*ldIter);
580 term = static_cast<RooArgSet*>(*lIter);
581 // PDF should be added to existing term if
582 // 1) It has overlapping normalization dependents with any other PDF in existing term
583 // 2) It has overlapping dependents of any class for which integration is requested
584 // 3) If normalization happens over multiple ranges, and those ranges are both defined
585 // in either observable
586
587 bool normOverlap = termNormDeps->overlaps(pdfNormDeps.begin(), pdfNormDeps.end());
588 //bool intOverlap = pdfIntSet->overlaps(*termAllDeps);
589
590 if (normOverlap) {
591// std::cout << GetName() << ": this term overlaps with term " << (*term) << " in normalization observables" << std::endl;
592
593 term->add(pdf);
594 termNormDeps->add(pdfNormDeps.begin(), pdfNormDeps.end(), false);
595 depAllList[j].add(pdfAllDeps.begin(), pdfAllDeps.end(), false);
596 if (termIntDeps) {
597 termIntDeps->add(pdfIntSet.begin(), pdfIntSet.end(), false);
598 }
599 if (termIntNoNormDeps) {
601 }
603 done = true;
604 }
605 }
606
607 // If not, create a new term
608 if (!done) {
609 if (!(pdfNormDeps.empty() && pdfAllDeps.empty() &&
610 pdfIntSet.empty()) || normSet.empty()) {
611 term = new RooArgSet("term");
612 RooArgSet *termNormDeps = new RooArgSet("termNormDeps");
613 depAllList.emplace_back(pdfAllDeps.begin(), pdfAllDeps.end(), "termAllDeps");
614 termIntDeps = new RooArgSet(pdfIntSet.begin(), pdfIntSet.end(), "termIntDeps");
615 depIntNoNormList.emplace_back(pdfIntNoNormDeps.begin(), pdfIntNoNormDeps.end(), "termIntNoNormDeps");
617
618 term->add(pdf);
619 termNormDeps->add(pdfNormDeps.begin(), pdfNormDeps.end(), false);
620
621 factorized.terms.Add(term);
622 factorized.norms.Add(termNormDeps);
623 factorized.ints.Add(termIntDeps);
624 }
625 }
626
627 }
628
629 // Loop over list of terms again to determine 'imported' observables
630 int i = 0;
632 auto lIter = factorized.terms.begin();
633 auto ldIter = factorized.norms.begin();
634 for(;lIter != factorized.terms.end(); (++lIter, ++ldIter, ++i)) {
635 normDeps = static_cast<RooArgSet*>(*ldIter);
636 term = static_cast<RooArgSet*>(*lIter);
637 // Make list of wholly imported dependents
639 impDeps.remove(*normDeps, true, true);
640 auto snap = new RooArgSet;
641 impDeps.snapshot(*snap);
642 factorized.imps.Add(snap);
643// std::cout << GetName() << ": list of imported dependents for term " << (*term) << " set to " << impDeps << std::endl ;
644
645 // Make list of cross dependents (term is self contained for these dependents,
646 // but components import dependents from other components)
647 auto crossDeps = std::unique_ptr<RooAbsCollection>{depIntNoNormList[i].selectCommon(*normDeps)};
648 snap = new RooArgSet;
649 crossDeps->snapshot(*snap);
650 factorized.cross.Add(snap);
651// std::cout << GetName() << ": list of cross dependents for term " << (*term) << " set to " << *crossDeps << std::endl ;
652 }
653}
654
655
656
657
658////////////////////////////////////////////////////////////////////////////////
659/// Return list of (partial) integrals of product terms for integration
660/// of p.d.f over observables iset while normalization over observables nset.
661/// Also return list of normalization sets to be used to evaluate
662/// each component in the list correctly.
663
665{
666 // Check if this configuration was created before
667 Int_t sterileIdx(-1);
668
669 if (static_cast<CacheElem*>(_cacheMgr.getObj(nset,iset,&sterileIdx,isetRangeName))) {
670 return _cacheMgr.lastIndex();
671 }
672
673 std::unique_ptr<CacheElem> cache = createCacheElem(nset, iset, isetRangeName);
674
675 // Store the partial integral list and return the assigned code
676 return _cacheMgr.setObj(nset, iset, cache.release(), RooNameReg::ptr(isetRangeName));
677}
678
679
680
681std::unique_ptr<RooProdPdf::CacheElem> RooProdPdf::createCacheElem(const RooArgSet* nset,
682 const RooArgSet* iset,
683 const char* isetRangeName) const
684{
685// std::cout << " FOLKERT::RooProdPdf::getPartIntList(" << GetName() <<") nset = " << (nset?*nset:RooArgSet()) << std::endl
686// << " _normRange = " << _normRange << std::endl
687// << " iset = " << (iset?*iset:RooArgSet()) << std::endl
688// << " isetRangeName = " << (isetRangeName?isetRangeName:"<null>") << std::endl ;
689
690 // Create containers for partial integral components to be generated
691 auto cache = std::make_unique<CacheElem>();
692
693 // Factorize the product in irreducible terms for this nset
695
696 // Normalization set used for factorization
697 RooArgSet factNset(nset ? (*nset) : _defNormSet);
698
700
701 // Group irriducible terms that need to be (partially) integrated together
702 std::list<std::vector<RooArgSet*>> groupedList;
705
706 // Loop over groups
707// std::cout<<"FK: pdf("<<GetName()<<") Starting selecting F(x|y)!"<< std::endl;
708 // Find groups of type F(x|y), i.e. termImpSet!=0, construct ratio object
709 std::map<std::string, RooArgSet> ratioTerms;
710 for (auto const& group : groupedList) {
711 if (1 == group.size()) {
712// std::cout<<"FK: Starting Single Term"<< std::endl;
713
714 RooArgSet* term = group[0];
715
716 Int_t termIdx = factorized.terms.IndexOf(term);
717 RooArgSet *norm=static_cast<RooArgSet*>(factorized.norms.At(termIdx));
718 RooArgSet *imps=static_cast<RooArgSet*>(factorized.imps.At(termIdx));
720 RooArgSet termImpSet(*imps);
721
722 // std::cout<<"FK: termImpSet.size() = "<<termImpSet.size()<< " " << termImpSet << std::endl;
723 // std::cout<<"FK: _refRangeName = "<<_refRangeName<< std::endl;
724
725 if (!termImpSet.empty() && nullptr != _refRangeName) {
726
727 // WVE we can skip this if the ref range is equal to the normalization range
728 // LM : avoid making integral ratio if range is the same. Why was not included ??? (same at line 857)
730// std::cout << "PREPARING RATIO HERE (SINGLE TERM)" << std::endl ;
732 std::ostringstream str; termImpSet.printValue(str);
733// std::cout << GetName() << "inserting ratio term" << std::endl;
734 ratioTerms[str.str()].addOwned(std::move(ratio));
735 }
736 }
737
738 } else {
739// std::cout<<"FK: Starting Composite Term"<< std::endl;
740
741 for (auto const& term : group) {
742
743 Int_t termIdx = factorized.terms.IndexOf(term);
744 RooArgSet *norm=static_cast<RooArgSet*>(factorized.norms.At(termIdx));
745 RooArgSet *imps=static_cast<RooArgSet*>(factorized.imps.At(termIdx));
747 RooArgSet termImpSet(*imps);
748
749 if (!termImpSet.empty() && nullptr != _refRangeName) {
750
751 // WVE we can skip this if the ref range is equal to the normalization range
753// std::cout << "PREPARING RATIO HERE (COMPOSITE TERM)" << std::endl ;
755 std::ostringstream str; termImpSet.printValue(str);
756 ratioTerms[str.str()].addOwned(std::move(ratio));
757 }
758 }
759 }
760 }
761
762 }
763
764 // Find groups with y as termNSet
765 // Replace G(y) with (G(y),ratio)
766 for (auto const& group : groupedList) {
767 for (auto const& term : group) {
768 Int_t termIdx = factorized.terms.IndexOf(term);
769 RooArgSet *norm = static_cast<RooArgSet*>(factorized.norms.At(termIdx));
770 RooArgSet *imps = static_cast<RooArgSet*>(factorized.imps.At(termIdx));
772 RooArgSet termImpSet(*imps);
773
774 // If termNset matches index of ratioTerms, insert ratio here
775 ostringstream str; termNSet.printValue(str);
776 if (!ratioTerms[str.str()].empty()) {
777// std::cout << "MUST INSERT RATIO OBJECT IN TERM (COMPOSITE)" << *term << std::endl;
778 term->add(ratioTerms[str.str()]);
779 cache->_ownedList.addOwned(std::move(ratioTerms[str.str()]));
780 }
781 }
782 }
783
784 for (auto const& group : groupedList) {
785// std::cout << GetName() << ":now processing group" << std::endl;
786// group->Print("1");
787
788 if (1 == group.size()) {
789// std::cout << "processing atomic item" << std::endl;
790 RooArgSet* term = group[0];
791
792 Int_t termIdx = factorized.terms.IndexOf(term);
793 RooArgSet *norm = factorized.termNormDeps(termIdx);
794 RooArgSet *integ = factorized.termIntDeps(termIdx);
795 RooArgSet *xdeps = factorized.termCrossDeps(termIdx);
796 RooArgSet *imps = factorized.termImpDeps(termIdx);
797
798 // Take list of normalization, integrated dependents, and
799 // cross-imported integrated dependents from factorization algorithm
803 RooArgSet termImpSet{*imps};
804
805 // Add prefab term to partIntList.
807 if (func.x0) {
808 cache->_partList.add(*func.x0);
809 if (func.isOwned) cache->_ownedList.addOwned(std::unique_ptr<RooAbsArg>{func.x0});
810
811 cache->_normList.emplace_back(std::make_unique<RooArgSet>());
812 norm->snapshot(*cache->_normList.back(), false);
813
814 cache->_numList.addOwned(std::move(func.x1));
815 cache->_denList.addOwned(std::move(func.x2));
816 }
817 } else {
818// std::cout << "processing composite item" << std::endl;
823 for (auto const &term : group) {
824 // std::cout << GetName() << ": processing term " << (*term) << " of composite item" << std::endl ;
825 Int_t termIdx = factorized.terms.IndexOf(term);
826 RooArgSet *norm = factorized.termNormDeps(termIdx);
827 RooArgSet *integ = factorized.termIntDeps(termIdx);
828 RooArgSet *xdeps = factorized.termCrossDeps(termIdx);
829 RooArgSet *imps = factorized.termImpDeps(termIdx);
830
834 RooArgSet termImpSet{*imps};
835
836 // Remove outer integration dependents from termISet
837 termISet.remove(outerIntDeps, true, true);
838
839 auto func = processProductTerm(nset, iset, isetRangeName, term, termNSet, termISet, true);
840 // std::cout << GetName() << ": created composite term component " << func.x0->GetName() << std::endl;
841 if (func.x0) {
842 compTermSet.add(*func.x0);
843 if (func.isOwned) cache->_ownedList.addOwned(std::unique_ptr<RooAbsArg>{func.x0});
844 compTermNorm.add(*norm, false);
845
846 compTermNum.add(*func.x1.release());
847 compTermDen.add(*func.x2.release());
848 //cache->_numList.add(*func.x1);
849 //cache->_denList.add(*func.x2);
850
851 }
852 }
853
854// std::cout << GetName() << ": constructing special composite product" << std::endl;
855// std::cout << GetName() << ": compTermSet = " ; compTermSet.Print("1");
856
857 // WVE THIS NEEDS TO BE REARRANGED
858
859 // compTermset is set van partial integrals to be multiplied
860 // prodtmp = product (compTermSet)
861 // inttmp = int ( prodtmp ) d (outerIntDeps) _range_isetRangeName
862
863 const std::string prodname = makeRGPPName("SPECPROD", compTermSet, outerIntDeps, RooArgSet(), isetRangeName);
864 auto prodtmp = std::make_unique<RooProduct>(prodname.c_str(), prodname.c_str(), compTermSet);
865
866 const std::string intname = makeRGPPName("SPECINT", compTermSet, outerIntDeps, RooArgSet(), isetRangeName);
867 auto inttmp = std::make_unique<RooRealIntegral>(intname.c_str(), intname.c_str(), *prodtmp, outerIntDeps, nullptr, nullptr, isetRangeName);
868 inttmp->setStringAttribute("PROD_TERM_TYPE", "SPECINT");
869
870 cache->_partList.add(*inttmp);
871
872 // Product of numerator terms
873 const string prodname_num = makeRGPPName("SPECPROD_NUM", compTermNum, RooArgSet(), RooArgSet(), nullptr);
874 auto prodtmp_num = std::make_unique<RooProduct>(prodname_num.c_str(), prodname_num.c_str(), compTermNum);
875 prodtmp_num->addOwnedComponents(compTermNum);
876
877 // Product of denominator terms
878 const string prodname_den = makeRGPPName("SPECPROD_DEN", compTermDen, RooArgSet(), RooArgSet(), nullptr);
879 auto prodtmp_den = std::make_unique<RooProduct>(prodname_den.c_str(), prodname_den.c_str(), compTermDen);
880 prodtmp_den->addOwnedComponents(compTermDen);
881
882 // Ratio
883 std::string name = Form("SPEC_RATIO(%s,%s)", prodname_num.c_str(), prodname_den.c_str());
884 auto ndr = std::make_unique<RooFormulaVar>(name.c_str(), "@0/@1", RooArgList(*prodtmp_num, *prodtmp_den));
885
886 // Integral of ratio
887 std::unique_ptr<RooAbsReal> numtmp{ndr->createIntegral(outerIntDeps,isetRangeName)};
888 numtmp->addOwnedComponents(std::move(ndr));
889
890 cache->_ownedList.addOwned(std::move(prodtmp));
891 cache->_ownedList.addOwned(std::move(inttmp));
892 cache->_ownedList.addOwned(std::move(prodtmp_num));
893 cache->_ownedList.addOwned(std::move(prodtmp_den));
894 cache->_numList.addOwned(std::move(numtmp));
895 cache->_denList.addOwned(std::unique_ptr<RooAbsArg>{static_cast<RooAbsArg*>(RooFit::RooConst(1).clone("1"))});
896 cache->_normList.emplace_back(std::make_unique<RooArgSet>());
897 compTermNorm.snapshot(*cache->_normList.back(), false);
898 }
899 }
900
901 // Need to rearrange product in case of multiple ranges
902 if (_normRange.Contains(",")) {
903 rearrangeProduct(*cache);
904 }
905
906 return cache;
907}
908
909
911{
912 // We own contents of all lists filled by factorizeProduct()
913 terms.Delete();
914 ints.Delete();
915 imps.Delete();
916 norms.Delete();
917 cross.Delete();
918}
919
920
921////////////////////////////////////////////////////////////////////////////////
922/// For single normalization ranges
923
924std::unique_ptr<RooAbsReal> RooProdPdf::makeCondPdfRatioCorr(RooAbsReal& pdf, const RooArgSet& termNset, const RooArgSet& /*termImpSet*/, const char* normRangeTmp, const char* refRange) const
925{
926 std::unique_ptr<RooAbsReal> ratio_num{pdf.createIntegral(termNset,normRangeTmp)};
927 std::unique_ptr<RooAbsReal> ratio_den{pdf.createIntegral(termNset,refRange)};
928 auto ratio = std::make_unique<RooFormulaVar>(Form("ratio(%s,%s)",ratio_num->GetName(),ratio_den->GetName()),"@0/@1",
930
931 ratio->addOwnedComponents(std::move(ratio_num));
932 ratio->addOwnedComponents(std::move(ratio_den));
933 ratio->setAttribute("RATIO_TERM") ;
934 return ratio ;
935}
936
937
938
939
940////////////////////////////////////////////////////////////////////////////////
941
943{
945
946 std::vector<std::string> rangeComps = ROOT::Split(_normRange.Data(), ",", /*skipEmpty=*/true);
947
948 std::map<std::string,RooArgSet> denListList ;
950 string specIntRange ;
951
952// std::cout << "THIS IS REARRANGEPRODUCT" << std::endl ;
953
954 for (std::size_t i = 0; i < cache._partList.size(); i++) {
955
956 RooAbsReal *part = static_cast<RooAbsReal*>(cache._partList.at(i));
957 RooAbsReal *num = static_cast<RooAbsReal*>(cache._numList.at(i));
958 RooAbsReal *den = static_cast<RooAbsReal*>(cache._denList.at(i));
959 i++;
960
961// std::cout << "now processing part " << part->GetName() << " of type " << part->getStringAttribute("PROD_TERM_TYPE") << std::endl ;
962// std::cout << "corresponding numerator = " << num->GetName() << std::endl ;
963// std::cout << "corresponding denominator = " << den->GetName() << std::endl ;
964
965
966 RooFormulaVar* ratio(nullptr) ;
968
969 if (string("SPECINT")==part->getStringAttribute("PROD_TERM_TYPE")) {
970
971 RooRealIntegral* orig = static_cast<RooRealIntegral*>(num);
972 auto specratio = static_cast<RooFormulaVar const*>(&orig->integrand()) ;
973 RooProduct* func = static_cast<RooProduct*>(specratio->getParameter(0)) ;
974
975 std::unique_ptr<RooArgSet> components{orig->getComponents()};
976 for(RooAbsArg * carg : *components) {
977 if (carg->getAttribute("RATIO_TERM")) {
978 ratio = static_cast<RooFormulaVar*>(carg) ;
979 break ;
980 }
981 }
982
983 if (ratio) {
984 RooCustomizer cust(*func,"blah") ;
985 cust.replaceArg(*ratio,RooFit::RooConst(1)) ;
986 nomList.add(*cust.build()) ;
987 } else {
988 nomList.add(*func) ;
989 }
990
991
992 } else {
993
994 // Find the ratio term
995 RooAbsReal* func = num;
996 // If top level object is integral, navigate to integrand
998 func = const_cast<RooAbsReal*>(&static_cast<RooRealIntegral*>(func)->integrand());
999 }
1000 if (func->InheritsFrom(RooProduct::Class())) {
1001// std::cout << "product term found: " ; func->Print() ;
1002 for(RooAbsArg * arg : static_cast<RooProduct*>(func)->components()) {
1003 if (arg->getAttribute("RATIO_TERM")) {
1004 ratio = static_cast<RooFormulaVar*>(arg) ;
1005 } else {
1006 origNumTerm.add(*arg) ;
1007 }
1008 }
1009 }
1010
1011 if (ratio) {
1012// std::cout << "Found ratio term in numerator: " << ratio->GetName() << std::endl ;
1013// std::cout << "Adding only original term to numerator: " << origNumTerm << std::endl ;
1014 nomList.add(origNumTerm) ;
1015 } else {
1016 nomList.add(*num) ;
1017 }
1018
1019 }
1020
1021 for (auto iter = rangeComps.begin() ; iter != rangeComps.end() ; ++iter) {
1022 // If denominator is an integral, make a clone with the integration range adjusted to
1023 // the selected component of the normalization integral
1024// std::cout << "NOW PROCESSING DENOMINATOR " << den->ClassName() << "::" << den->GetName() << std::endl ;
1025
1026 if (string("SPECINT")==part->getStringAttribute("PROD_TERM_TYPE")) {
1027
1028// std::cout << "create integral: SPECINT case" << std::endl ;
1029 RooRealIntegral* orig = static_cast<RooRealIntegral*>(num);
1030 auto specRatio = static_cast<RooFormulaVar const*>(&orig->integrand()) ;
1031 specIntDeps.add(orig->intVars()) ;
1032 if (orig->intRange()) {
1033 specIntRange = orig->intRange() ;
1034 }
1035 //RooProduct* numtmp = (RooProduct*) specRatio->getParameter(0) ;
1036 RooProduct* dentmp = static_cast<RooProduct*>(specRatio->getParameter(1)) ;
1037
1038// std::cout << "numtmp = " << numtmp->ClassName() << "::" << numtmp->GetName() << std::endl ;
1039// std::cout << "dentmp = " << dentmp->ClassName() << "::" << dentmp->GetName() << std::endl ;
1040
1041// std::cout << "denominator components are " << dentmp->components() << std::endl ;
1042 for (auto* parg : static_range_cast<RooAbsReal*>(dentmp->components())) {
1043// std::cout << "now processing denominator component " << parg->ClassName() << "::" << parg->GetName() << std::endl ;
1044
1045 if (ratio && parg->dependsOn(*ratio)) {
1046// std::cout << "depends in value of ratio" << std::endl ;
1047
1048 // Make specialize ratio instance
1049 std::unique_ptr<RooAbsReal> specializedRatio{specializeRatio(*(RooFormulaVar*)ratio,iter->c_str())};
1050
1051// std::cout << "specRatio = " << std::endl ;
1052// specializedRatio->printComponentTree() ;
1053
1054 // Replace generic ratio with specialized ratio
1055 RooAbsArg *partCust(nullptr) ;
1056 if (parg->InheritsFrom(RooAddition::Class())) {
1057
1058
1059
1060 RooAddition* tmpadd = static_cast<RooAddition*>(parg) ;
1061
1062 RooCustomizer cust(*tmpadd->list1().first(),Form("blah_%s",iter->c_str())) ;
1063 cust.replaceArg(*ratio,*specializedRatio) ;
1064 partCust = cust.build() ;
1065
1066 } else {
1067 RooCustomizer cust(*parg,Form("blah_%s",iter->c_str())) ;
1068 cust.replaceArg(*ratio,*specializedRatio) ;
1069 partCust = cust.build() ;
1070 }
1071
1072 // Print customized denominator
1073// std::cout << "customized function = " << std::endl ;
1074// partCust->printComponentTree() ;
1075
1076 std::unique_ptr<RooAbsReal> specializedPartCust{specializeIntegral(*static_cast<RooAbsReal*>(partCust),iter->c_str())};
1077
1078 // Finally divide again by ratio
1079 string name = Form("%s_divided_by_ratio",specializedPartCust->GetName()) ;
1080 auto specIntFinal = std::make_unique<RooFormulaVar>(name.c_str(),"@0/@1",RooArgList(*specializedPartCust,*specializedRatio)) ;
1081 specIntFinal->addOwnedComponents(std::move(specializedPartCust));
1082 specIntFinal->addOwnedComponents(std::move(specializedRatio));
1083
1084 denListList[*iter].addOwned(std::move(specIntFinal));
1085 } else {
1086
1087// std::cout << "does NOT depend on value of ratio" << std::endl ;
1088// parg->Print("t") ;
1089
1090 denListList[*iter].addOwned(specializeIntegral(*parg,iter->c_str()));
1091
1092 }
1093 }
1094// std::cout << "end iteration over denominator components" << std::endl ;
1095 } else {
1096
1097 if (ratio) {
1098
1099 std::unique_ptr<RooAbsReal> specRatio{specializeRatio(*(RooFormulaVar*)ratio,iter->c_str())};
1100
1101 // If integral is 'Int r(y)*g(y) dy ' then divide a posteriori by r(y)
1102// std::cout << "have ratio, orig den = " << den->GetName() << std::endl ;
1103
1104 RooArgSet tmp(origNumTerm) ;
1105 tmp.add(*specRatio) ;
1106 const string pname = makeRGPPName("PROD",tmp,RooArgSet(),RooArgSet(),nullptr) ;
1107 auto specDenProd = std::make_unique<RooProduct>(pname.c_str(),pname.c_str(),tmp) ;
1108 std::unique_ptr<RooAbsReal> specInt;
1109
1110 if (den->InheritsFrom(RooRealIntegral::Class())) {
1111 specInt = std::unique_ptr<RooAbsReal>{specDenProd->createIntegral((static_cast<RooRealIntegral*>(den))->intVars(),iter->c_str())};
1112 specInt->addOwnedComponents(std::move(specDenProd));
1113 } else if (den->InheritsFrom(RooAddition::Class())) {
1114 RooAddition* orig = static_cast<RooAddition*>(den) ;
1115 RooRealIntegral* origInt = static_cast<RooRealIntegral*>(orig->list1().first()) ;
1116 specInt = std::unique_ptr<RooAbsReal>{specDenProd->createIntegral(origInt->intVars(),iter->c_str())};
1117 specInt->addOwnedComponents(std::move(specDenProd));
1118 } else {
1119 throw string("this should not happen") ;
1120 }
1121
1122 //RooAbsReal* specInt = specializeIntegral(*den,iter->c_str()) ;
1123 string name = Form("%s_divided_by_ratio",specInt->GetName()) ;
1124 auto specIntFinal = std::make_unique<RooFormulaVar>(name.c_str(),"@0/@1",RooArgList(*specInt,*specRatio)) ;
1125 specIntFinal->addOwnedComponents(std::move(specInt));
1126 specIntFinal->addOwnedComponents(std::move(specRatio));
1127 denListList[*iter].addOwned(std::move(specIntFinal));
1128 } else {
1129 denListList[*iter].addOwned(specializeIntegral(*den,iter->c_str()));
1130 }
1131
1132 }
1133 }
1134
1135 }
1136
1137 // Do not rearrange terms if numerator and denominator are effectively empty
1138 if (nomList.empty()) {
1139 return ;
1140 }
1141
1142 string name = Form("%s_numerator",GetName()) ;
1143 // WVE FIX THIS (2)
1144
1145 std::unique_ptr<RooAbsReal> numerator = std::make_unique<RooProduct>(name.c_str(),name.c_str(),nomList) ;
1146
1148// std::cout << "nomList = " << nomList << std::endl ;
1149 for (map<string,RooArgSet>::iterator iter = denListList.begin() ; iter != denListList.end() ; ++iter) {
1150// std::cout << "denList[" << iter->first << "] = " << iter->second << std::endl ;
1151 name = Form("%s_denominator_comp_%s",GetName(),iter->first.c_str()) ;
1152 // WVE FIX THIS (2)
1153 RooProduct* prod_comp = new RooProduct(name.c_str(),name.c_str(),iter->second) ;
1154 prod_comp->addOwnedComponents(std::move(iter->second));
1155 products.add(*prod_comp) ;
1156 }
1157 name = Form("%s_denominator_sum",GetName()) ;
1158 RooAbsReal* norm = new RooAddition(name.c_str(),name.c_str(),products) ;
1159 norm->addOwnedComponents(products) ;
1160
1161 if (!specIntDeps.empty()) {
1162 // Apply posterior integration required for SPECINT case
1163
1164 string namesr = Form("SPEC_RATIO(%s,%s)",numerator->GetName(),norm->GetName()) ;
1165 RooFormulaVar* ndr = new RooFormulaVar(namesr.c_str(),"@0/@1",RooArgList(*numerator,*norm)) ;
1166 ndr->addOwnedComponents(std::move(numerator));
1167
1168 // Integral of ratio
1169 numerator = std::unique_ptr<RooAbsReal>{ndr->createIntegral(specIntDeps,specIntRange.c_str())};
1170
1171 norm = static_cast<RooAbsReal*>(RooFit::RooConst(1).Clone()) ;
1172 }
1173
1174
1175// std::cout << "numerator" << std::endl ;
1176// numerator->printComponentTree("",0,5) ;
1177// std::cout << "denominator" << std::endl ;
1178// norm->printComponentTree("",0,5) ;
1179
1180
1181 // WVE DEBUG
1182 //RooMsgService::instance().debugWorkspace()->import(RooArgSet(*numerator,*norm)) ;
1183
1184 cache._rearrangedNum = std::move(numerator);
1185 cache._rearrangedDen.reset(norm);
1186 cache._isRearranged = true ;
1187
1188}
1189
1190
1191////////////////////////////////////////////////////////////////////////////////
1192
1193std::unique_ptr<RooAbsReal> RooProdPdf::specializeRatio(RooFormulaVar& input, const char* targetRangeName) const
1194{
1195 RooRealIntegral* numint = static_cast<RooRealIntegral*>(input.getParameter(0)) ;
1196 RooRealIntegral* denint = static_cast<RooRealIntegral*>(input.getParameter(1)) ;
1197
1198 std::unique_ptr<RooAbsReal> numint_spec{specializeIntegral(*numint,targetRangeName)};
1199
1200 std::unique_ptr<RooAbsReal> ret = std::make_unique<RooFormulaVar>(Form("ratio(%s,%s)",numint_spec->GetName(),denint->GetName()),"@0/@1",RooArgList(*numint_spec,*denint)) ;
1201 ret->addOwnedComponents(std::move(numint_spec));
1202
1203 return ret;
1204}
1205
1206
1207
1208////////////////////////////////////////////////////////////////////////////////
1209
1210std::unique_ptr<RooAbsReal> RooProdPdf::specializeIntegral(RooAbsReal& input, const char* targetRangeName) const
1211{
1212 if (input.InheritsFrom(RooRealIntegral::Class())) {
1213
1214 // If input is integral, recreate integral but override integration range to be targetRangeName
1215 RooRealIntegral* orig = static_cast<RooRealIntegral*>(&input) ;
1216// std::cout << "creating integral: integrand = " << orig->integrand().GetName() << " vars = " << orig->intVars() << " range = " << targetRangeName << std::endl ;
1217 return std::unique_ptr<RooAbsReal>{orig->integrand().createIntegral(orig->intVars(),targetRangeName)};
1218
1219 } else if (input.InheritsFrom(RooAddition::Class())) {
1220
1221 // If input is sum of integrals, recreate integral from first component of set, but override integration range to be targetRangeName
1222 RooAddition* orig = static_cast<RooAddition*>(&input) ;
1223 RooRealIntegral* origInt = static_cast<RooRealIntegral*>(orig->list1().first()) ;
1224// std::cout << "creating integral from addition: integrand = " << origInt->integrand().GetName() << " vars = " << origInt->intVars() << " range = " << targetRangeName << std::endl ;
1225 return std::unique_ptr<RooAbsReal>{origInt->integrand().createIntegral(origInt->intVars(),targetRangeName)};
1226 }
1227
1228 std::stringstream errMsg;
1229 errMsg << "specializeIntegral: unknown input type " << input.ClassName() << "::" << input.GetName();
1230 throw std::runtime_error(errMsg.str());
1231}
1232
1233
1234////////////////////////////////////////////////////////////////////////////////
1235/// Group product into terms that can be calculated independently
1236
1237void RooProdPdf::groupProductTerms(std::list<std::vector<RooArgSet*>>& groupedTerms, RooArgSet& outerIntDeps,
1238 Factorized const &factorized) const
1239{
1240 // Start out with each term in its own group
1241 for(auto * term : static_range_cast<RooArgSet*>(factorized.terms)) {
1242 groupedTerms.emplace_back();
1243 groupedTerms.back().emplace_back(term) ;
1244 }
1245
1246 // Make list of imported dependents that occur in any term
1249 allImpDeps.add(*impDeps,false) ;
1250 }
1251
1252 // Make list of integrated dependents that occur in any term
1255 allIntDeps.add(*intDeps,false) ;
1256 }
1257
1258 outerIntDeps.removeAll() ;
1259 outerIntDeps.add(*std::unique_ptr<RooArgSet>{allIntDeps.selectCommon(allImpDeps)});
1260
1261 // Now iteratively merge groups that should be (partially) integrated together
1263
1264 // Collect groups that feature this dependent
1265 std::vector<RooArgSet*>* newGroup = nullptr ;
1266
1267 // Loop over groups
1268 bool needMerge = false ;
1269 auto group = groupedTerms.begin();
1270 auto nGroups = groupedTerms.size();
1271 for (size_t iGroup = 0; iGroup < nGroups; ++iGroup) {
1272
1273 // See if any term in this group depends in any ay on outerDepInt
1274 for (auto const& term2 : *group) {
1275
1276 Int_t termIdx = factorized.terms.IndexOf(term2) ;
1277 if (factorized.termNormDeps(termIdx)->contains(*outerIntDep) ||
1278 factorized.termIntDeps(termIdx)->contains(*outerIntDep) ||
1279 factorized.termImpDeps(termIdx)->contains(*outerIntDep)) {
1280 needMerge = true ;
1281 }
1282
1283 }
1284
1285 if (needMerge) {
1286 // Create composite group if not yet existing
1287 if (newGroup==nullptr) {
1288 groupedTerms.emplace_back() ;
1289 newGroup = &groupedTerms.back() ;
1290 }
1291
1292 // Add terms of this group to new term
1293 for (auto& term2 : *group) {
1294 newGroup->emplace_back(term2) ;
1295 }
1296
1297 // Remove this non-owning group from list
1298 group = groupedTerms.erase(group);
1299 } else {
1300 ++group;
1301 }
1302 }
1303
1304 }
1305}
1306
1307
1308
1309////////////////////////////////////////////////////////////////////////////////
1310/// Calculate integrals of factorized product terms over observables iset while normalized
1311/// to observables in nset.
1312
1314 const RooArgSet* term,const RooArgSet& termNSet, const RooArgSet& termISet,
1315 bool forceWrap) const
1316{
1318
1319 // CASE I: factorizing term: term is integrated over all normalizing observables
1320 // -----------------------------------------------------------------------------
1321 // Check if all observbales of this term are integrated. If so the term cancels
1322 if (!termNSet.empty() && termNSet.size()==termISet.size() && isetRangeName==nullptr) {
1323 // Term factorizes
1324 return ret ;
1325 }
1326
1327 // CASE II: Dropped terms: if term is entirely unnormalized, it should be dropped
1328 // ------------------------------------------------------------------------------
1329 if (nset && termNSet.empty()) {
1330 // Drop terms that are not asked to be normalized
1331 return ret ;
1332 }
1333
1334 if (iset && !termISet.empty()) {
1335 if (term->size()==1) {
1336
1337 // CASE IIIa: Normalized and partially integrated single PDF term
1338 //---------------------------------------------------------------
1339
1340 RooAbsPdf* pdf = static_cast<RooAbsPdf*>(term->first()) ;
1341
1342 ret.x0 = std::unique_ptr<RooAbsReal>{pdf->createIntegral(termISet,termNSet,isetRangeName)}.release();
1343 ret.x0->setOperMode(operMode()) ;
1344 ret.x0->setStringAttribute("PROD_TERM_TYPE","IIIa") ;
1345
1346 ret.isOwned=true ;
1347
1348 // Split mode results
1349 ret.x1 = std::unique_ptr<RooAbsReal>{pdf->createIntegral(termISet,isetRangeName)};
1350 ret.x2 = std::unique_ptr<RooAbsReal>{pdf->createIntegral(termNSet,normRange())};
1351
1352 return ret ;
1353
1354 } else {
1355
1356 // CASE IIIb: Normalized and partially integrated composite PDF term
1357 //---------------------------------------------------------------
1358
1359 // Use auxiliary class RooGenProdProj to calculate this term
1360 const std::string name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
1361 ret.x0 = new RooGenProdProj(name.c_str(),name.c_str(),*term,termISet,termNSet,isetRangeName) ;
1362 ret.x0->setStringAttribute("PROD_TERM_TYPE","IIIb") ;
1363 ret.x0->setOperMode(operMode()) ;
1364
1365 ret.isOwned=true ;
1366
1367 const std::string name1 = makeRGPPName("PROD",*term,RooArgSet(),RooArgSet(),nullptr) ;
1368
1369 // WVE FIX THIS
1370 RooProduct* tmp_prod = new RooProduct(name1.c_str(),name1.c_str(),*term) ;
1371
1372 ret.x1 = std::unique_ptr<RooAbsReal>{tmp_prod->createIntegral(termISet,isetRangeName)};
1373 ret.x2 = std::unique_ptr<RooAbsReal>{tmp_prod->createIntegral(termNSet,normRange())};
1374
1375 return ret ;
1376 }
1377 }
1378
1379 // CASE IVa: Normalized non-integrated composite PDF term
1380 // -------------------------------------------------------
1381 if (nset && !nset->empty() && term->size()>1) {
1382 // Composite term needs normalized integration
1383
1384 const std::string name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
1385 ret.x0 = new RooGenProdProj(name.c_str(),name.c_str(),*term,termISet,termNSet,isetRangeName,normRange()) ;
1386 ret.x0->setExpensiveObjectCache(expensiveObjectCache()) ;
1387
1388 ret.x0->setStringAttribute("PROD_TERM_TYPE","IVa") ;
1389 ret.x0->setOperMode(operMode()) ;
1390
1391 ret.isOwned=true ;
1392
1393 const std::string name1 = makeRGPPName("PROD",*term,RooArgSet(),RooArgSet(),nullptr) ;
1394
1395 // WVE FIX THIS
1396 RooProduct* tmp_prod = new RooProduct(name1.c_str(),name1.c_str(),*term) ;
1397
1398 ret.x1 = std::unique_ptr<RooAbsReal>{tmp_prod->createIntegral(termISet,isetRangeName)};
1399 ret.x2 = std::unique_ptr<RooAbsReal>{tmp_prod->createIntegral(termNSet,normRange())};
1400
1401 return ret ;
1402 }
1403
1404 // CASE IVb: Normalized, non-integrated single PDF term
1405 // -----------------------------------------------------
1406 for (auto* pdf : static_range_cast<RooAbsPdf*>(*term)) {
1407
1408 ret.isOwned = false;
1409 RooAbsReal *ret0 = pdf;
1410
1411 if (forceWrap) {
1412 ret.isOwned = true;
1413 // Construct representative name of normalization wrapper
1414 std::string name = pdf->GetName() + ("_NORM[" + RooHelpers::getColonSeparatedNameString(termNSet, ','));
1415 name += normRange() ? ('|' + std::string{normRange()} + ']') : "]";
1416 ret0 = new RooRealIntegral(name.c_str(),name.c_str(),*pdf,RooArgSet(),&termNSet);
1417 }
1418
1419 ret0->setStringAttribute("PROD_TERM_TYPE","IVb") ;
1420 ret.x0 = ret0;
1421 ret.x1 = std::unique_ptr<RooAbsReal>{pdf->createIntegral(RooArgSet())};
1422 ret.x2 = std::unique_ptr<RooAbsReal>{pdf->createIntegral(termNSet,normRange())};
1423 return ret;
1424 }
1425
1426 coutE(Eval) << "RooProdPdf::processProductTerm(" << GetName() << ") unidentified term!!!" << std::endl ;
1427 return ret ;
1428}
1429
1430
1431
1432
1433////////////////////////////////////////////////////////////////////////////////
1434/// Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
1435
1436std::string RooProdPdf::makeRGPPName(const char* pfx, const RooArgSet& term, const RooArgSet& iset,
1437 const RooArgSet& nset, const char* isetRangeName) const
1438{
1439 // Make an appropriate automatic name for a RooGenProdProj object in getPartIntList()
1440
1441 std::ostringstream os(pfx);
1442 os << "[";
1443
1444 // Encode component names
1445 bool first(true) ;
1446 for (auto const* pdf : static_range_cast<RooAbsPdf*>(term)) {
1447 if (!first) os << "_X_";
1448 first = false;
1449 os << pdf->GetName();
1450 }
1451 os << "]" << integralNameSuffix(iset,&nset,isetRangeName,true);
1452
1453 return os.str();
1454}
1455
1456
1457
1458////////////////////////////////////////////////////////////////////////////////
1459/// Force RooRealIntegral to offer all observables for internal integration
1460
1462{
1463 return true ;
1464}
1465
1466
1467
1468////////////////////////////////////////////////////////////////////////////////
1469/// Determine which part (if any) of given integral can be performed analytically.
1470/// If any analytical integration is possible, return integration scenario code.
1471///
1472/// RooProdPdf implements two strategies in implementing analytical integrals
1473///
1474/// First, PDF components whose entire set of dependents are requested to be integrated
1475/// can be dropped from the product, as they will integrate out to 1 by construction
1476///
1477/// Second, RooProdPdf queries each remaining component PDF for its analytical integration
1478/// capability of the requested set ('allVars'). It finds the largest common set of variables
1479/// that can be integrated by all remaining components. If such a set exists, it reconfirms that
1480/// each component is capable of analytically integrating the common set, and combines the components
1481/// individual integration codes into a single integration code valid for RooProdPdf.
1482
1484 const RooArgSet* normSet, const char* rangeName) const
1485{
1486 if (_forceNumInt) return 0 ;
1487
1488 // Declare that we can analytically integrate all requested observables
1489 analVars.add(allVars) ;
1490
1491 // Retrieve (or create) the required partial integral list
1492 Int_t code = getPartIntList(normSet,&allVars,rangeName);
1493
1494 return code+1 ;
1495}
1496
1497
1498
1499
1500////////////////////////////////////////////////////////////////////////////////
1501/// Return analytical integral defined by given scenario code
1502
1503double RooProdPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
1504{
1505 // No integration scenario
1506 if (code==0) {
1507 return getVal(normSet) ;
1508 }
1509
1510
1511 // WVE needs adaptation for rangename feature
1512
1513 // Partial integration scenarios
1514 CacheElem* cache = static_cast<CacheElem*>(_cacheMgr.getObjByIndex(code-1)) ;
1515
1516 // If cache has been sterilized, revive this slot
1517 if (cache==nullptr) {
1518 std::unique_ptr<RooArgSet> vars{getParameters(RooArgSet())} ;
1519 RooArgSet nset = _cacheMgr.selectFromSet1(*vars, code-1) ;
1520 RooArgSet iset = _cacheMgr.selectFromSet2(*vars, code-1) ;
1521
1522 Int_t code2 = getPartIntList(&nset, &iset, rangeName) ;
1523
1524 // preceding call to getPartIntList guarantees non-null return
1525 // coverity[NULL_RETURNS]
1526 cache = static_cast<CacheElem*>(_cacheMgr.getObj(&nset,&iset,&code2,rangeName)) ;
1527 }
1528
1529 double val = calculate(*cache,true) ;
1530// std::cout << "RPP::aIWN(" << GetName() << ") ,code = " << code << ", value = " << val << std::endl ;
1531
1532 return val ;
1533}
1534
1535
1536
1537////////////////////////////////////////////////////////////////////////////////
1538/// If this product contains exactly one extendable p.d.f return the extension abilities of
1539/// that p.d.f, otherwise return CanNotBeExtended
1540
1542{
1543 return (_extendedIndex>=0) ? (static_cast<RooAbsPdf*>(_pdfList.at(_extendedIndex)))->extendMode() : CanNotBeExtended ;
1544}
1545
1546
1547
1548////////////////////////////////////////////////////////////////////////////////
1549/// Return the expected number of events associated with the extendable input PDF
1550/// in the product. If there is no extendable term, abort.
1551
1552double RooProdPdf::expectedEvents(const RooArgSet* nset) const
1553{
1554 if (_extendedIndex<0) {
1555 coutF(Generation) << "Requesting expected number of events from a RooProdPdf that does not contain an extended p.d.f" << std::endl ;
1556 throw std::logic_error(std::string("RooProdPdf ") + GetName() + " could not be extended.");
1557 }
1558
1559 return static_cast<RooAbsPdf*>(_pdfList.at(_extendedIndex))->expectedEvents(nset) ;
1560}
1561
1562std::unique_ptr<RooAbsReal> RooProdPdf::createExpectedEventsFunc(const RooArgSet* nset) const
1563{
1564 if (_extendedIndex<0) {
1565 coutF(Generation) << "Requesting expected number of events from a RooProdPdf that does not contain an extended p.d.f" << std::endl ;
1566 throw std::logic_error(std::string("RooProdPdf ") + GetName() + " could not be extended.");
1567 }
1568
1569 return static_cast<RooAbsPdf*>(_pdfList.at(_extendedIndex))->createExpectedEventsFunc(nset);
1570}
1571
1572
1573////////////////////////////////////////////////////////////////////////////////
1574/// Return generator context optimized for generating events from product p.d.f.s
1575
1577 const RooArgSet* auxProto, bool verbose) const
1578{
1579 if (_useDefaultGen) return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ;
1580 return new RooProdGenContext(*this,vars,prototype,auxProto,verbose) ;
1581}
1582
1583
1584
1585////////////////////////////////////////////////////////////////////////////////
1586/// Query internal generation capabilities of component p.d.f.s and aggregate capabilities
1587/// into master configuration passed to the generator context
1588
1590{
1591 if (!_useDefaultGen) return 0 ;
1592
1593 // Find the subset directVars that only depend on a single PDF in the product
1595 for (auto const* arg : directVars) {
1596 if (isDirectGenSafe(*arg)) directSafe.add(*arg) ;
1597 }
1598
1599
1600 // Now find direct integrator for relevant components ;
1601 std::vector<Int_t> code;
1602 code.reserve(64);
1603 for (auto const* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
1605 Int_t pdfCode = pdf->getGenerator(directSafe,pdfDirect,staticInitOK);
1606 code.push_back(pdfCode);
1607 if (pdfCode != 0) {
1608 generateVars.add(pdfDirect) ;
1609 }
1610 }
1611
1612
1613 if (!generateVars.empty()) {
1614 Int_t masterCode = _genCode.store(code) ;
1615 return masterCode+1 ;
1616 } else {
1617 return 0 ;
1618 }
1619}
1620
1621
1622
1623////////////////////////////////////////////////////////////////////////////////
1624/// Forward one-time initialization call to component generation initialization
1625/// methods.
1626
1628{
1629 if (!_useDefaultGen) return ;
1630
1631 const std::vector<Int_t>& codeList = _genCode.retrieve(code-1) ;
1632 Int_t i(0) ;
1633 for (auto* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
1634 if (codeList[i]!=0) {
1635 pdf->initGenerator(codeList[i]) ;
1636 }
1637 i++ ;
1638 }
1639}
1640
1641
1642
1643////////////////////////////////////////////////////////////////////////////////
1644/// Generate a single event with configuration specified by 'code'
1645/// Defer internal generation to components as encoded in the _genCode
1646/// registry for given generator code.
1647
1649{
1650 if (!_useDefaultGen) return ;
1651
1652 const std::vector<Int_t>& codeList = _genCode.retrieve(code-1) ;
1653 Int_t i(0) ;
1654 for (auto* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
1655 if (codeList[i]!=0) {
1656 pdf->generateEvent(codeList[i]) ;
1657 }
1658 i++ ;
1659 }
1660}
1661
1662
1663
1664////////////////////////////////////////////////////////////////////////////////
1665/// Return RooAbsArg components contained in the cache
1666
1668{
1669 RooArgList ret ;
1670 ret.add(_partList) ;
1671 ret.add(_numList) ;
1672 ret.add(_denList) ;
1673 if (_rearrangedNum) ret.add(*_rearrangedNum) ;
1674 if (_rearrangedDen) ret.add(*_rearrangedDen) ;
1675 return ret ;
1676
1677}
1678
1679
1680
1681////////////////////////////////////////////////////////////////////////////////
1682/// Hook function to print cache contents in tree printing of RooProdPdf
1683
1685{
1686 if (curElem==0) {
1687 os << indent << "RooProdPdf begin partial integral cache" << std::endl ;
1688 }
1689
1690 auto indent2 = std::string(indent) + "[" + std::to_string(curElem) + "]";
1691 for(auto const& arg : _partList) {
1692 arg->printCompactTree(os,indent2.c_str()) ;
1693 }
1694
1695 if (curElem==maxElem) {
1696 os << indent << "RooProdPdf end partial integral cache" << std::endl ;
1697 }
1698}
1699
1700
1701
1702////////////////////////////////////////////////////////////////////////////////
1703/// Forward determination of safety of internal generator code to
1704/// component p.d.f that would generate the given observable
1705
1707{
1708 // Only override base class behaviour if default generator method is enabled
1709 if (!_useDefaultGen) return RooAbsPdf::isDirectGenSafe(arg) ;
1710
1711 // Argument may appear in only one PDF component
1712 RooAbsPdf* thePdf(nullptr) ;
1713 for (auto* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
1714
1715 if (pdf->dependsOn(arg)) {
1716 // Found PDF depending on arg
1717
1718 // If multiple PDFs depend on arg directGen is not safe
1719 if (thePdf) return false ;
1720
1721 thePdf = pdf ;
1722 }
1723 }
1724 // Forward call to relevant component PDF
1725 return thePdf?(thePdf->isDirectGenSafe(arg)):false ;
1726}
1727
1728
1729
1730////////////////////////////////////////////////////////////////////////////////
1731/// Look up user specified normalization set for given input PDF component
1732
1734{
1735 Int_t idx = _pdfList.index(&pdf) ;
1736 if (idx<0) return nullptr;
1737 return _pdfNSetList[idx].get() ;
1738}
1739
1740
1741
1742/// Add some full PDFs to the factors of this RooProdPdf.
1744{
1745 size_t numExtended = (_extendedIndex==-1) ? 0 : 1;
1746
1747 for(auto arg : pdfs) {
1748 RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg);
1749 if (!pdf) {
1750 coutW(InputArguments) << "RooProdPdf::addPdfs(" << GetName() << ") list arg "
1751 << arg->GetName() << " is not a PDF, ignored" << std::endl ;
1752 continue;
1753 }
1754 if(pdf->canBeExtended()) {
1755 if (_extendedIndex == -1) {
1757 } else {
1758 numExtended++;
1759 }
1760 }
1761 _pdfList.add(*pdf);
1762 _pdfNSetList.emplace_back(std::make_unique<RooArgSet>("nset"));
1763 }
1764
1765 // Protect against multiple extended terms
1766 if (numExtended>1) {
1767 coutW(InputArguments) << "RooProdPdf::addPdfs(" << GetName()
1768 << ") WARNING: multiple components with extended terms detected,"
1769 << " product will not be extendable." << std::endl ;
1770 _extendedIndex = -1 ;
1771 }
1772
1773 // Reset cache
1774 _cacheMgr.reset() ;
1775
1776}
1777
1778/// Remove some PDFs from the factors of this RooProdPdf.
1780{
1781 // Remember what the extended PDF is
1782 RooAbsArg const* extPdf = _extendedIndex >= 0 ? &_pdfList[_extendedIndex] : nullptr;
1783
1784 // Actually remove the PDFs and associated nsets
1785 for(size_t i=0;i < _pdfList.size(); i++) {
1786 if(pdfs.contains(_pdfList[i])) {
1788 _pdfNSetList.erase(_pdfNSetList.begin()+i);
1789 i--;
1790 }
1791 }
1792
1793 // Since we may have removed PDFs from the list, the index of the extended
1794 // PDF in the list needs to be updated. The new index might also be -1 if the
1795 // extended PDF got removed.
1796 if(extPdf) {
1798 }
1799
1800 // Reset cache
1801 _cacheMgr.reset() ;
1802}
1803
1804
1805namespace {
1806
1807std::vector<TNamed const*> sortedNamePtrs(RooAbsCollection const& col)
1808{
1809 std::vector<TNamed const*> ptrs;
1810 ptrs.reserve(col.size());
1811 for(RooAbsArg* arg : col) {
1812 ptrs.push_back(arg->namePtr());
1813 }
1814 std::sort(ptrs.begin(), ptrs.end());
1815 return ptrs;
1816}
1817
1818bool sortedNamePtrsOverlap(std::vector<TNamed const*> const& ptrsA, std::vector<TNamed const*> const& ptrsB)
1819{
1820 auto pA = ptrsA.begin();
1821 auto pB = ptrsB.begin();
1822 while (pA != ptrsA.end() && pB != ptrsB.end()) {
1823 if (*pA < *pB) {
1824 ++pA;
1825 } else if (*pB < *pA) {
1826 ++pB;
1827 } else {
1828 return true;
1829 }
1830 }
1831 return false;
1832}
1833
1834} // namespace
1835
1836
1837////////////////////////////////////////////////////////////////////////////////
1838/// Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
1839/// The observables set is required to distinguish unambiguously p.d.f in terms
1840/// of observables and parameters, which are not constraints, and p.d.fs in terms
1841/// of parameters only, which can serve as constraints p.d.f.s
1842/// The pdfParams output parameter communicates to the caller which parameter
1843/// are used in the pdfs that are not constraints.
1844
1846{
1847 auto constraints = new RooArgSet{"constraints"};
1848
1849 // For the optimized implementation of checking if two collections overlap by name.
1850 auto observablesNamePtrs = sortedNamePtrs(observables);
1852
1853 // Loop over PDF components
1854 for (std::size_t iPdf = 0; iPdf < _pdfList.size(); ++iPdf) {
1855 auto * pdf = static_cast<RooAbsPdf*>(&_pdfList[iPdf]);
1856
1857 RooArgSet tmp;
1858 pdf->getParameters(nullptr, tmp);
1859
1860 // A constraint term is a p.d.f that doesn't contribute to the
1861 // expectedEvents() and does not depend on any of the listed observables
1862 // but does depends on any of the parameters that should be constrained
1863 bool isConstraint = false;
1864
1865 if(static_cast<int>(iPdf) != _extendedIndex) {
1866 auto tmpNamePtrs = sortedNamePtrs(tmp);
1867 // Before, there were calls to `pdf->dependsOn()` here, but they were very
1868 // expensive for large computation graphs! Given that we have to traverse
1869 // the computation graph with a call to `pdf->getParameters()` anyway, we
1870 // can just check if the set of all variables operlaps with the observables
1871 // or constraind parameters.
1872 //
1873 // We are using an optimized implementation of overlap checking. Because
1874 // the overlap is checked by name, we can check overlap of the
1875 // corresponding name pointers. The optimization can't be in
1876 // RooAbsCollection itself, because it is crucial that the memory for the
1877 // non-tmp name pointers is not reallocated for each pdf.
1880 }
1881 if (isConstraint) {
1882 constraints->add(*pdf) ;
1883 } else {
1884 // We only want to add parameter, not observables. Since a call like
1885 // `pdf->getParameters(&observables)` would be expensive, we take the set
1886 // of all variables and remove the ovservables, which is much cheaper. In
1887 // a call to `pdf->getParameters(&observables)`, the observables are
1888 // matched by name, so we have to pass the `matchByNameOnly` here.
1889 tmp.remove(observables, /*silent=*/false, /*matchByNameOnly=*/true);
1890 pdfParams.add(tmp,true) ;
1891 }
1892 }
1893
1894 return constraints;
1895}
1896
1897
1898
1899
1900////////////////////////////////////////////////////////////////////////////////
1901/// Return all parameter constraint p.d.f.s on parameters listed in constrainedParams.
1902/// The observables set is required to distinguish unambiguously p.d.f in terms
1903/// of observables and parameters, which are not constraints, and p.d.fs in terms
1904/// of parameters only, which can serve as constraints p.d.f.s
1905
1907{
1908 RooArgSet* connectedPars = new RooArgSet("connectedPars") ;
1909 for (std::size_t iPdf = 0; iPdf < _pdfList.size(); ++iPdf) {
1910 auto * pdf = static_cast<RooAbsPdf*>(&_pdfList[iPdf]);
1911 // Check if term is relevant, either because it provides a propablity
1912 // density in the observables or because it is used for the expected
1913 // events.
1914 if (static_cast<int>(iPdf) == _extendedIndex || pdf->dependsOn(observables)) {
1915 RooArgSet tmp;
1916 pdf->getParameters(&observables, tmp);
1917 connectedPars->add(tmp) ;
1918 }
1919 }
1920 return connectedPars ;
1921}
1922
1923
1924
1925
1926////////////////////////////////////////////////////////////////////////////////
1927
1929{
1931 if (!nset || nset->empty()) return ;
1932
1933 // Get/create appropriate term list for this normalization set
1934 Int_t code = getPartIntList(nset, nullptr);
1935 RooArgList & plist = static_cast<CacheElem*>(_cacheMgr.getObj(nset, &code))->_partList;
1936
1937 // Strip any terms from params that do not depend on any term
1939 for (auto param : *params) {
1940 bool anyDep(false) ;
1941 for (auto term : plist) {
1942 if (term->dependsOnValue(*param)) {
1943 anyDep=true ;
1944 }
1945 }
1946 if (!anyDep) {
1947 tostrip.add(*param) ;
1948 }
1949 }
1950
1951 if (!tostrip.empty()) {
1952 params->remove(tostrip,true,true);
1953 }
1954
1955}
1956
1957
1958
1959////////////////////////////////////////////////////////////////////////////////
1960/// Interface function used by test statistics to freeze choice of range
1961/// for interpretation of conditional product terms
1962
1963void RooProdPdf::selectNormalizationRange(const char* rangeName, bool force)
1964{
1965 if (!force && _refRangeName) {
1966 return ;
1967 }
1968
1969 fixRefRange(rangeName) ;
1970}
1971
1972
1973
1974
1975////////////////////////////////////////////////////////////////////////////////
1976
1977void RooProdPdf::fixRefRange(const char* rangeName)
1978{
1979 _refRangeName = const_cast<TNamed*>(RooNameReg::ptr(rangeName));
1980}
1981
1982
1983
1984////////////////////////////////////////////////////////////////////////////////
1985/// Forward the plot sampling hint from the p.d.f. that defines the observable obs
1986
1987std::list<double>* RooProdPdf::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
1988{
1989 for (auto const* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
1990 if (std::list<double>* hint = pdf->plotSamplingHint(obs,xlo,xhi)) {
1991 return hint ;
1992 }
1993 }
1994
1995 return nullptr;
1996}
1997
1998
1999
2000////////////////////////////////////////////////////////////////////////////////
2001/// If all components that depend on obs are binned that so is the product
2002
2004{
2005 for (auto const* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
2006 if (pdf->dependsOn(obs) && !pdf->isBinnedDistribution(obs)) {
2007 return false ;
2008 }
2009 }
2010
2011 return true ;
2012}
2013
2014
2015
2016
2017
2018
2019////////////////////////////////////////////////////////////////////////////////
2020/// Forward the plot sampling hint from the p.d.f. that defines the observable obs
2021
2022std::list<double>* RooProdPdf::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
2023{
2024 for (auto const* pdf : static_range_cast<RooAbsPdf*>(_pdfList)) {
2025 if (std::list<double>* hint = pdf->binBoundaries(obs,xlo,xhi)) {
2026 return hint ;
2027 }
2028 }
2029
2030 return nullptr;
2031}
2032
2033
2034////////////////////////////////////////////////////////////////////////////////
2035/// Label OK'ed components of a RooProdPdf with cache-and-track, _and_ label all RooProdPdf
2036/// descendants with extra information about (conditional) normalization, needed to be able
2037/// to Cache-And-Track them outside the RooProdPdf context.
2038
2040{
2041 for (const auto parg : _pdfList) {
2042
2043 if (parg->canNodeBeCached()==Always) {
2044 trackNodes.add(*parg) ;
2045// std::cout << "tracking node RooProdPdf component " << parg << " " << parg->ClassName() << "::" << parg->GetName() << std::endl ;
2046
2047 // Additional processing to fix normalization sets in case product defines conditional observables
2048 if (RooArgSet* pdf_nset = findPdfNSet(static_cast<RooAbsPdf&>(*parg))) {
2049 // Check if conditional normalization is specified
2051 if (string("nset")==pdf_nset->GetName() && !pdf_nset->empty()) {
2052 parg->setStringAttribute("CATNormSet",getColonSeparatedNameString(*pdf_nset).c_str()) ;
2053 }
2054 if (string("cset")==pdf_nset->GetName()) {
2055 parg->setStringAttribute("CATCondSet",getColonSeparatedNameString(*pdf_nset).c_str()) ;
2056 }
2057 } else {
2058 coutW(Optimization) << "RooProdPdf::setCacheAndTrackHints(" << GetName() << ") WARNING product pdf does not specify a normalization set for component " << parg->GetName() << std::endl ;
2059 }
2060 }
2061 }
2062}
2063
2064
2065
2066////////////////////////////////////////////////////////////////////////////////
2067/// Customized printing of arguments of a RooProdPdf to more intuitively reflect the contents of the
2068/// product operator construction
2069
2070void RooProdPdf::printMetaArgs(ostream& os) const
2071{
2072 for (std::size_t i=0 ; i<_pdfList.size() ; i++) {
2073 if (i>0) os << " * " ;
2074 RooArgSet* ncset = _pdfNSetList[i].get() ;
2075 os << _pdfList.at(i)->GetName() ;
2076 if (!ncset->empty()) {
2077 if (string("nset")==ncset->GetName()) {
2078 os << *ncset ;
2079 } else {
2080 os << "|" ;
2081 bool first(true) ;
2082 for (auto const* arg : *ncset) {
2083 if (!first) {
2084 os << "," ;
2085 } else {
2086 first = false ;
2087 }
2088 os << arg->GetName() ;
2089 }
2090 }
2091 }
2092 }
2093 os << " " ;
2094}
2095
2096
2097
2098////////////////////////////////////////////////////////////////////////////////
2099/// Implement support for node removal
2100
2102{
2103 if (nameChange && _pdfList.find("REMOVAL_DUMMY")) {
2104
2105 cxcoutD(LinkStateMgmt) << "RooProdPdf::redirectServersHook(" << GetName() << "): removing REMOVAL_DUMMY" << std::endl ;
2106
2107 // Remove node from _pdfList proxy and remove corresponding entry from normset list
2108 RooAbsArg* pdfDel = _pdfList.find("REMOVAL_DUMMY") ;
2109
2110 _pdfNSetList.erase(_pdfNSetList.begin() + _pdfList.index("REMOVAL_DUMMY")) ;
2112
2113 // Clear caches
2114 _cacheMgr.reset() ;
2115 }
2116
2117 // If the replaced server is an observable that is used in any of the
2118 // normalization sets for conditional fits, replace the element in the
2119 // normalization set too.
2120 for(std::unique_ptr<RooArgSet> const& normSet : _pdfNSetList) {
2121 for(RooAbsArg * arg : *normSet) {
2122 if(RooAbsArg * newArg = arg->findNewServer(newServerList, nameChange)) {
2123 // Since normSet is owning, the original arg is now deleted.
2124 normSet->replace(arg, std::unique_ptr<RooAbsArg>{newArg->cloneTree()});
2125 }
2126 }
2127 }
2128
2130}
2131
2132void RooProdPdf::CacheElem::writeToStream(std::ostream& os) const {
2133 using namespace RooHelpers;
2134 os << "_partList\n";
2135 os << getColonSeparatedNameString(_partList) << "\n";
2136 os << "_numList\n";
2137 os << getColonSeparatedNameString(_numList) << "\n";
2138 os << "_denList\n";
2139 os << getColonSeparatedNameString(_denList) << "\n";
2140 os << "_ownedList\n";
2141 os << getColonSeparatedNameString(_ownedList) << "\n";
2142 os << "_normList\n";
2143 for(auto const& set : _normList) {
2144 os << getColonSeparatedNameString(*set) << "\n";
2145 }
2146 os << "_isRearranged" << "\n";
2147 os << _isRearranged << "\n";
2148 os << "_rearrangedNum" << "\n";
2149 if(_rearrangedNum) {
2150 os << getColonSeparatedNameString(*_rearrangedNum) << "\n";
2151 } else {
2152 os << "nullptr" << "\n";
2153 }
2154 os << "_rearrangedDen" << "\n";
2155 if(_rearrangedDen) {
2156 os << getColonSeparatedNameString(*_rearrangedDen) << "\n";
2157 } else {
2158 os << "nullptr" << "\n";
2159 }
2160}
2161
2162std::unique_ptr<RooArgSet> RooProdPdf::fillNormSetForServer(RooArgSet const &normSet, RooAbsArg const &server) const
2163{
2164 if (normSet.empty())
2165 return nullptr;
2166 auto *pdfNset = findPdfNSet(static_cast<RooAbsPdf const &>(server));
2167 if (pdfNset && !pdfNset->empty()) {
2168 std::unique_ptr<RooArgSet> out;
2169 if (0 == strcmp("cset", pdfNset->GetName())) {
2170 // If the name of the normalization set is "cset", it doesn't contain the
2171 // normalization set but the conditional observables that should *not* be
2172 // normalized over.
2173 out = std::make_unique<RooArgSet>(normSet);
2175 out->selectCommon(*pdfNset, common);
2176 out->remove(common);
2177 } else {
2178 out = std::make_unique<RooArgSet>(*pdfNset);
2179 }
2180 // prefix also the arguments in the normSets if they have not already been
2181 if (auto prefix = getStringAttribute("__prefix__")) {
2182 for (RooAbsArg *arg : *out) {
2183 if (!arg->getStringAttribute("__prefix__")) {
2184 arg->SetName((std::string(prefix) + arg->GetName()).c_str());
2185 arg->setStringAttribute("__prefix__", prefix);
2186 }
2187 }
2188 }
2189 return out;
2190 } else {
2191 return nullptr;
2192 }
2193}
2194
2195std::unique_ptr<RooAbsArg>
2197{
2198 if (ctx.likelihoodMode()) {
2199 auto binnedInfo = RooHelpers::getBinnedL(*this);
2200 if (binnedInfo.binnedPdf && binnedInfo.binnedPdf != this) {
2201 return binnedInfo.binnedPdf->compileForNormSet(normSet, ctx);
2202 }
2203 }
2204
2205 std::unique_ptr<RooProdPdf> prodPdfClone{static_cast<RooProdPdf *>(this->Clone())};
2207
2208 for (const auto server : prodPdfClone->servers()) {
2210 RooArgSet const &nset = nsetForServer ? *nsetForServer : normSet;
2211
2213 server->getObservables(&nset, depList);
2214
2216 }
2217
2218 auto fixedProdPdf = std::make_unique<RooFit::Detail::RooFixedProdPdf>(std::move(prodPdfClone), normSet);
2220
2221 return fixedProdPdf;
2222}
2223
2224namespace RooFit {
2225namespace Detail {
2226
2227RooFixedProdPdf::RooFixedProdPdf(std::unique_ptr<RooProdPdf> &&prodPdf, RooArgSet const &normSet)
2228 : RooAbsPdf(prodPdf->GetName(), prodPdf->GetTitle()),
2229 _normSet{normSet},
2230 _servers("!servers", "List of servers", this),
2231 _prodPdf{std::move(prodPdf)}
2232{
2233 initialize();
2234}
2235
2237 : RooAbsPdf(other, name),
2238 _normSet{other._normSet},
2239 _servers("!servers", "List of servers", this),
2240 _prodPdf{static_cast<RooProdPdf *>(other._prodPdf->Clone())}
2241{
2242 initialize();
2243}
2244
2246{
2247 _cache = _prodPdf->createCacheElem(&_normSet, nullptr);
2248 auto &cache = *_cache;
2249
2250 // The actual servers for a given normalization set depend on whether the
2251 // cache is rearranged or not. See RooProdPdf::calculateBatch to see
2252 // which args in the cache are used directly.
2253 if (cache._isRearranged) {
2256 } else {
2257 for (std::size_t i = 0; i < cache._partList.size(); ++i) {
2259 }
2260 }
2261}
2262
2263} // namespace Detail
2264} // namespace RooFit
ROOT::RRangeCast< T, false, Range_t > static_range_cast(Range_t &&coll)
#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...
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) 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:220
RooFixedProdPdf(std::unique_ptr< RooProdPdf > &&prodPdf, RooArgSet const &normSet)
std::unique_ptr< RooProdPdf > _prodPdf
Definition RooProdPdf.h:274
RooProdPdf::CacheElem const & cache() const
Definition RooProdPdf.h:264
std::unique_ptr< RooProdPdf::CacheElem > _cache
Definition RooProdPdf.h:272
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...
General form of projected integral of product of PDFs, utility class for RooProdPdf.
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:196
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:199
std::vector< std::unique_ptr< RooArgSet > > _pdfNSetList
List of PDF component normalization sets.
Definition RooProdPdf.h:195
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:201
RooAICRegistry _genCode
! Registry of composite direct generator codes
Definition RooProdPdf.h:191
void addPdfs(RooAbsCollection const &pdfs)
Add some full PDFs to the factors of this RooProdPdf.
RooListProxy _pdfList
List of PDF components.
Definition RooProdPdf.h:194
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:185
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:204
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 doEvalImpl(RooAbsArg const *caller, const RooProdPdf::CacheElem &cache, RooFit::EvalContext &) const
Evaluate product of PDFs in batch mode.
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:193
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:47
Mother of all ROOT objects.
Definition TObject.h:41
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition TObject.cxx:542
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={})
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:64
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