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