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