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