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