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