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