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