/*****************************************************************************
* Project: RooFit *
* Package: RooFitCore *
* File: $Id: RooProdPdf.cc,v 1.60 2005/06/23 07:37:30 wverkerke Exp $
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
// -- CLASS DESCRIPTION [PDF] --
// RooProdPdf is an efficient implementation of a product of PDFs of the form
//
// PDF_1 * PDF_2 * ... * PDF_N
//
// PDFs may share observables. If that is the case any irreducable subset
// of PDFS that share observables will be normalized with explicit numeric
// integration as any built-in normalization will no longer be valid.
//
// Alternatively, products using conditional PDFs can be defined, e.g.
//
// F(x|y) * G(y)
//
// meaning a pdf F(x) _given_ y and a PDF G(y). In this contruction F is only
// normalized w.r.t x and G is normalized w.r.t y. The product in this construction
// is properly normalized.
//
// If exactly one of the component PDFs supports extended likelihood fits, the
// product will also be usable in extended mode, returning the number of expected
// events from the extendable component PDF. The extendable component does not
// have to appear in any specific place in the list.
//
#include "RooFit.h"
#include "TIterator.h"
#include "TIterator.h"
#include "RooProdPdf.h"
#include "RooRealProxy.h"
#include "RooProdGenContext.h"
#include "RooGenProdProj.h"
#include "RooProduct.h"
#include "RooNameReg.h"
ClassImp(RooProdPdf)
;
RooProdPdf::RooProdPdf(const char *name, const char *title, Double_t cutOff) :
RooAbsPdf(name,title),
_partListMgr(10),
_partOwnedListMgr(10),
_genCode(10),
_cutOff(cutOff),
_pdfList("_pdfList","List of PDFs",this),
_extendedIndex(-1),
_useDefaultGen(kFALSE)
{
// Dummy constructor
}
RooProdPdf::RooProdPdf(const char *name, const char *title,
RooAbsPdf& pdf1, RooAbsPdf& pdf2, Double_t cutOff) :
RooAbsPdf(name,title),
_partListMgr(10),
_partOwnedListMgr(10),
_genCode(10),
_cutOff(cutOff),
_pdfList("_pdfList","List of PDFs",this),
_pdfIter(_pdfList.createIterator()),
_extendedIndex(-1),
_useDefaultGen(kFALSE)
{
// Constructor with 2 PDFs (most frequent use case).
//
// The optional cutOff parameter can be used as a speed optimization if
// one or more of the PDF have sizable regions with very small values,
// which would pull the entire product of PDFs to zero in those regions.
//
// After each PDF multiplication, the running product is compared with
// the cutOff parameter. If the running product is smaller than the
// cutOff value, the product series is terminated and remaining PDFs
// are not evaluated.
//
// There is no magic value of the cutOff, the user should experiment
// to find the appropriate balance between speed and precision.
// If a cutoff is specified, the PDFs most likely to be small should
// be put first in the product. The default cutOff value is zero.
_pdfList.add(pdf1) ;
RooArgSet* nset1 = new RooArgSet("nset1") ;
_pdfNSetList.Add(nset1) ;
if (pdf1.canBeExtended()) {
_extendedIndex = _pdfList.index(&pdf1) ;
}
_pdfList.add(pdf2) ;
RooArgSet* nset2 = new RooArgSet("nset2") ;
_pdfNSetList.Add(nset2) ;
if (pdf2.canBeExtended()) {
if (_extendedIndex>=0) {
// Protect against multiple extended terms
cout << "RooProdPdf::RooProdPdf(" << GetName()
<< ") multiple components with extended terms detected,"
<< " product will not be extendible." << endl ;
_extendedIndex=-1 ;
} else {
_extendedIndex=_pdfList.index(&pdf2) ;
}
}
}
RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgList& pdfList, Double_t cutOff) :
RooAbsPdf(name,title),
_partListMgr(10),
_partOwnedListMgr(10),
_genCode(10),
_cutOff(cutOff),
_pdfList("_pdfList","List of PDFs",this),
_pdfIter(_pdfList.createIterator()),
_extendedIndex(-1),
_useDefaultGen(kFALSE)
{
// Constructor from a list of PDFs
//
// The optional cutOff parameter can be used as a speed optimization if
// one or more of the PDF have sizable regions with very small values,
// which would pull the entire product of PDFs to zero in those regions.
//
// After each PDF multiplication, the running product is compared with
// the cutOff parameter. If the running product is smaller than the
// cutOff value, the product series is terminated and remaining PDFs
// are not evaluated.
//
// There is no magic value of the cutOff, the user should experiment
// to find the appropriate balance between speed and precision.
// If a cutoff is specified, the PDFs most likely to be small should
// be put first in the product. The default cutOff value is zero.
TIterator* iter = pdfList.createIterator() ;
RooAbsArg* arg ;
Int_t numExtended(0) ;
while((arg=(RooAbsArg*)iter->Next())) {
RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
if (!pdf) {
cout << "RooProdPdf::RooProdPdf(" << GetName() << ") list arg "
<< pdf->GetName() << " is not a PDF, ignored" << endl ;
continue ;
}
_pdfList.add(*pdf) ;
RooArgSet* nset = new RooArgSet("nset") ;
_pdfNSetList.Add(nset) ;
if (pdf->canBeExtended()) {
_extendedIndex = _pdfList.index(pdf) ;
numExtended++ ;
}
}
// Protect against multiple extended terms
if (numExtended>1) {
cout << "RooProdPdf::RooProdPdf(" << GetName()
<< ") WARNING: multiple components with extended terms detected,"
<< " product will not be extendible." << endl ;
_extendedIndex = -1 ;
}
delete iter ;
}
RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet,
const RooCmdArg& arg1, const RooCmdArg& arg2,
const RooCmdArg& arg3, const RooCmdArg& arg4,
const RooCmdArg& arg5, const RooCmdArg& arg6,
const RooCmdArg& arg7, const RooCmdArg& arg8) :
RooAbsPdf(name,title),
_partListMgr(10),
_partOwnedListMgr(10),
_genCode(10),
_cutOff(0),
_pdfList("_pdfList","List of PDFs",this),
_pdfIter(_pdfList.createIterator()),
_extendedIndex(-1),
_useDefaultGen(kFALSE)
{
// Constructor from named argument list
//
// fullPdf -- Set of 'regular' PDFS that are normalized over all their observables
// ConditionalPdf(pdfSet,depSet) -- Add PDF to product with condition that it
// only be normalized over specified observables
// any remaining observables will be conditional
// observables
//
//
// For example, given a PDF F(x,y) and G(y)
//
// RooProdPdf("P","P",G,Partial(F,x)) will construct a 2-dimensional PDF as follows:
//
// P(x,y) = G(y)/Int[y]G(y) * F(x,y)/Int[x]G(x,y)
//
// which is a well normalized and properly defined PDF, but different from the
//
// P'(x,y) = F(x,y)*G(y) / Int[x,y] F(x,y)*G(y)
//
// In the former case the y distribution of P is identical to that of G, while
// F only is used to determine the correlation between X and Y. In the latter
// case the Y distribution is defined by the product of F and G.
//
// This P(x,y) construction is analoguous to generating events from F(x,y) with
// a prototype dataset sampled from G(y)
RooLinkedList l ;
l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
initializeFromCmdArgList(fullPdfSet,l) ;
}
RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet, const RooLinkedList& cmdArgList) :
RooAbsPdf(name,title),
_partListMgr(10),
_partOwnedListMgr(10),
_genCode(10),
_cutOff(0),
_pdfList("_pdfList","List of PDFs",this),
_pdfIter(_pdfList.createIterator()),
_extendedIndex(-1),
_useDefaultGen(kFALSE)
{
// cout << "RooProdPdf::ctor" << endl ;
// cmdArgList.Print("v") ;
initializeFromCmdArgList(fullPdfSet, cmdArgList) ;
}
RooProdPdf::RooProdPdf(const RooProdPdf& other, const char* name) :
RooAbsPdf(other,name),
_partListMgr(other._partListMgr,kTRUE),
_partOwnedListMgr(other._partOwnedListMgr,kTRUE),
_genCode(other._genCode),
_cutOff(other._cutOff),
_pdfList("_pdfList",this,other._pdfList),
_pdfIter(_pdfList.createIterator()),
_extendedIndex(other._extendedIndex),
_useDefaultGen(other._useDefaultGen)
{
// Copy constructor
// Clone contents of normalizarion set list
TIterator* iter = other._pdfNSetList.MakeIterator() ;
RooArgSet* nset ;
while((nset=(RooArgSet*)iter->Next())) {
_pdfNSetList.Add(nset->snapshot()) ;
}
delete iter ;
}
void RooProdPdf::initializeFromCmdArgList(const RooArgSet& fullPdfSet, const RooLinkedList& l)
{
// Initialize RooProdPdf from a list of RooCmdArg configuration arguments
// Process set of full PDFS
TIterator* siter = fullPdfSet.createIterator() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)siter->Next())) {
_pdfList.add(*pdf) ;
RooArgSet* nset1 = new RooArgSet("nset1") ;
_pdfNSetList.Add(nset1) ;
}
delete siter ;
// Process list of conditional PDFs
TIterator* iter = l.MakeIterator() ;
RooCmdArg* carg ;
while((carg=(RooCmdArg*)iter->Next())) {
if (!TString(carg->GetName()).CompareTo("Conditional")) {
RooArgSet* pdfSet = (RooArgSet*) carg->getObject(0) ;
RooArgSet* normSet = (RooArgSet*) carg->getObject(1) ;
TIterator* siter = pdfSet->createIterator() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)siter->Next())) {
_pdfList.add(*pdf) ;
_pdfNSetList.Add(normSet->snapshot()) ;
}
delete siter ;
} else if (TString(carg->GetName()).CompareTo("")) {
cout << "Unknown arg: " << carg->GetName() << endl ;
}
}
delete iter ;
}
RooProdPdf::~RooProdPdf()
{
// Destructor
Int_t i ;
_pdfNSetList.Delete() ;
for (i=0 ; i<10 ; i++) {
_partNormListCache[i].Delete() ;
}
delete _pdfIter ;
}
Double_t RooProdPdf::getVal(const RooArgSet* set) const
{
_curNormSet = (RooArgSet*)set ;
return RooAbsPdf::getVal(set) ;
}
Double_t RooProdPdf::evaluate() const
{
// Calculate current unnormalized value of object
// WVE NSET may not exists, but getPartIntList may need it...
Int_t code ;
RooArgList *plist ;
RooLinkedList *nlist ;
getPartIntList(_curNormSet,0,plist,nlist,code) ;
return calculate(plist,nlist) ;
}
Double_t RooProdPdf::calculate(const RooArgList* partIntList, const RooLinkedList* normSetList) const
{
// Calculate running product of pdfs, skipping factorized components
RooAbsReal* partInt ;
RooArgSet* normSet ;
Double_t value(1.0) ;
Int_t n = partIntList->getSize() ;
Int_t i ;
for (i=0 ; i<n ; i++) {
partInt = ((RooAbsReal*)partIntList->at(i)) ;
normSet = ((RooArgSet*)normSetList->At(i)) ;
Double_t piVal = partInt->getVal(normSet->getSize()>0 ? normSet : 0) ;
value *= piVal ;
// if (_verboseEval<0) {
// cout << "RPP:calc(" << GetName() << "): value *= " << piVal << " (" << partInt->GetName() << ") nset = " ; normSet->Print("1") ;
// }
if (value<_cutOff) {
//cout << "RooProdPdf::calculate(" << GetName() << ") calculation cut off after " << partInt->GetName() << endl ;
break ;
}
}
return value ;
}
void RooProdPdf::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
// cout << "RooProdPdf::factorizeProduct(" << GetName() << ")" << endl ;
// cout << " normSet = " ; normSet.Print("1") ;
// cout << " intSet = " ; intSet.Print("1") ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
// List of all term dependents: normalization and imported
RooLinkedList depAllList ;
RooLinkedList depIntNoNormList ;
// Setup lists for factorization terms and their dependents
RooArgSet* term(0) ;
RooArgSet* termNormDeps(0) ;
RooArgSet* termAllDeps(0) ;
RooArgSet* termIntDeps(0) ;
RooArgSet* termIntNoNormDeps(0) ;
TIterator* lIter = termList.MakeIterator() ;
TIterator* ldIter = normList.MakeIterator() ;
TIterator* laIter = depAllList.MakeIterator() ;
TIterator* nIter = _pdfNSetList.MakeIterator() ;
RooArgSet* pdfNSet ;
// Loop over the PDFs
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
pdfNSet = (RooArgSet*) nIter->Next() ;
lIter->Reset() ;
ldIter->Reset() ;
laIter->Reset() ;
// Reduce pdfNSet to actual dependents
pdfNSet = pdf->getObservables(*pdfNSet) ;
// cout << "factorize: pdf = " << pdf->GetName() << " nset = " ; pdfNSet->Print("1") ;
RooArgSet pdfNormDeps ; // Dependents to be normalized for the PDF
RooArgSet pdfAllDeps ; // All dependents of this PDF ;
// Make list of all dependents of this PDF
RooArgSet* tmp = pdf->getObservables(normSet) ;
pdfAllDeps.add(*tmp) ;
delete tmp ;
// Make list of normalization dependents for this PDF ;
if (pdfNSet->getSize()>0) {
RooArgSet* tmp = (RooArgSet*) pdfAllDeps.selectCommon(*pdfNSet) ;
pdfNormDeps.add(*tmp) ;
delete tmp ;
} else {
pdfNormDeps.add(pdfAllDeps) ;
}
RooArgSet* pdfIntSet = pdf->getObservables(intSet) ;
RooArgSet pdfIntNoNormDeps(*pdfIntSet) ;
pdfIntNoNormDeps.remove(pdfNormDeps,kTRUE,kTRUE) ;
// cout << "factorize: all pdf dependents = " ; pdfAllDeps.Print("1") ;
// cout << "factorize: normalized pdf dependents = " ; pdfNormDeps.Print("1") ;
// cout << "factorize: integrated pdf dependents = " ; pdfIntSet->Print("1") ;
// cout << "factorize: int/no-norm pdf dependents = " ; pdfIntNoNormDeps.Print("1") ;
// Check if this PDF has dependents overlapping with one of the existing terms
Bool_t done(kFALSE) ;
while((term=(RooArgSet*)lIter->Next())) {
termNormDeps=(RooArgSet*)ldIter->Next() ;
termAllDeps=(RooArgSet*)laIter->Next() ;
// PDF should be added to existing term if
// 1) It has overlapping normalization dependents with any other PDF in existing term
// 2) It has overlapping dependents of any class for which integration is requested
Bool_t normOverlap = pdfNormDeps.overlaps(*termNormDeps) ;
// Bool_t intOverlap = pdfIntSet->overlaps(*termAllDeps) ;
//if (normOverlap || intOverlap) {
if (normOverlap) {
// cout << "adding pdf " << pdf->GetName() << " to term " ; term->Print("1") ;
// cout << " with termNormDeps " ; termNormDeps->Print("1") ;
term->add(*pdf) ;
termNormDeps->add(pdfNormDeps,kFALSE) ;
termAllDeps->add(pdfAllDeps,kFALSE) ;
termIntDeps->add(*pdfIntSet,kFALSE) ;
termIntNoNormDeps->add(pdfIntNoNormDeps,kFALSE) ;
done = kTRUE ;
}
}
// If not, create a new term
if (!done) {
// cout << "creating new term" << endl ;
term = new RooArgSet("term") ;
termNormDeps = new RooArgSet("termNormDeps") ;
termAllDeps = new RooArgSet("termAllDeps") ;
termIntDeps = new RooArgSet("termIntDeps") ;
termIntNoNormDeps = new RooArgSet("termIntNoNormDeps") ;
term->add(*pdf) ;
termNormDeps->add(pdfNormDeps,kFALSE) ;
termAllDeps->add(pdfAllDeps,kFALSE) ;
termIntDeps->add(*pdfIntSet,kFALSE) ;
termIntNoNormDeps->add(pdfIntNoNormDeps,kFALSE) ;
termList.Add(term) ;
normList.Add(termNormDeps) ;
depAllList.Add(termAllDeps) ;
intList.Add(termIntDeps) ;
depIntNoNormList.Add(termIntNoNormDeps) ;
}
// We own the reduced version of pdfNSet
delete pdfNSet ;
delete pdfIntSet ;
}
// Loop over list of terms again to determine 'imported' observables
lIter->Reset() ;
ldIter->Reset() ;
laIter->Reset() ;
TIterator* innIter = depIntNoNormList.MakeIterator() ;
while((term=(RooArgSet*)lIter->Next())) {
RooArgSet* normDeps = (RooArgSet*) ldIter->Next() ;
RooArgSet* allDeps = (RooArgSet*) laIter->Next() ;
RooArgSet* intNoNormDeps = (RooArgSet*) innIter->Next() ;
// Make list of wholly imported dependents
RooArgSet impDeps(*allDeps) ;
impDeps.remove(*normDeps,kTRUE,kTRUE) ;
impDepList.Add(impDeps.snapshot()) ;
// Make list of cross dependents (term is self contained for these dependents,
// but components import dependents from other components)
RooArgSet* crossDeps = (RooArgSet*) intNoNormDeps->selectCommon(*normDeps) ;
crossDepList.Add(crossDeps->snapshot()) ;
delete crossDeps ;
}
// lIter->Reset() ;
// cout << "list of terms:" << endl ;
// while((term=(RooArgSet*)lIter->Next())) {
// term->Print("1") ;
// }
depAllList.Delete() ;
depIntNoNormList.Delete() ;
delete nIter ;
delete lIter ;
delete ldIter ;
delete laIter ;
delete innIter ;
return ;
}
void RooProdPdf::getPartIntList(const RooArgSet* nset, const RooArgSet* iset,
pRooArgList& partList, pRooLinkedList& nsetList,
Int_t& code, const char* isetRangeName) const
{
// Check if this configuration was created before
Int_t sterileIdx(-1) ;
RooArgList* partIntList = _partListMgr.getNormList(this,nset,iset,&sterileIdx,RooNameReg::ptr(isetRangeName)) ;
if (partIntList) {
code = _partListMgr.lastIndex() ;
partList = partIntList ;
nsetList = _partNormListCache+code ;
return ;
}
// Create containers for partial integral components to be generated
partIntList = new RooArgList("partIntList") ;
RooArgList* partIntOwnedList = new RooArgList("partIntList") ;
RooLinkedList* partIntNormList = new RooLinkedList ;
// Factorize the product in irreducible terms for this nset
RooLinkedList terms, norms, imp, ints, cross ;
factorizeProduct(nset?(*nset):RooArgSet(),iset?(*iset):RooArgSet(),terms,norms,imp,cross,ints) ;
RooArgSet *norm, *integ, *xdeps ;
// Group irriducible terms that need to be (partially) integrated together
RooLinkedList groupedList ;
RooArgSet outerIntDeps ;
groupProductTerms(groupedList,outerIntDeps,terms,norms,imp,ints,cross) ;
TIterator* gIter = groupedList.MakeIterator() ;
RooLinkedList* group ;
while((group=(RooLinkedList*)gIter->Next())) {
// group->Print("1") ;
if (group->GetSize()==1) {
RooArgSet* term = (RooArgSet*) group->At(0) ;
Int_t termIdx = terms.IndexOf(term) ;
norm=(RooArgSet*) norms.At(termIdx) ;
integ=(RooArgSet*) ints.At(termIdx) ;
xdeps=(RooArgSet*)cross.At(termIdx) ;
RooArgSet termNSet, termISet, termXSet ;
// Take list of normalization, integrated dependents from factorization algorithm
termISet.add(*integ) ;
termNSet.add(*norm) ;
// Cross-imported integrated dependents
termXSet.add(*xdeps) ;
// Add prefab term to partIntList.
Bool_t isOwned ;
RooAbsReal* func = processProductTerm(nset,iset,isetRangeName,term,termNSet,termISet,isOwned) ;
if (func) {
partIntList->add(*func) ;
if (isOwned) partIntOwnedList->addOwned(*func) ;
partIntNormList->Add(norm->snapshot()) ;
}
} else {
// cout << "prod: composite group encountered!" << endl ;
RooArgSet compTermSet, compTermNorm ;
TIterator* tIter = group->MakeIterator() ;
RooArgSet* term ;
while((term=(RooArgSet*)tIter->Next())) {
Int_t termIdx = terms.IndexOf(term) ;
norm=(RooArgSet*) norms.At(termIdx) ;
integ=(RooArgSet*) ints.At(termIdx) ;
xdeps=(RooArgSet*)cross.At(termIdx) ;
RooArgSet termNSet, termISet, termXSet ;
termISet.add(*integ) ;
termNSet.add(*norm) ;
termXSet.add(*xdeps) ;
// Remove outer integration dependents from termISet
termISet.remove(outerIntDeps,kTRUE,kTRUE) ;
// cout << "termISet = " ; termISet.Print("1") ;
Bool_t isOwned ;
RooAbsReal* func = processProductTerm(nset,iset,isetRangeName,term,termNSet,termISet,isOwned,kTRUE) ;
// cout << "created composite term component " << func->GetName() << endl ;
if (func) {
compTermSet.add(*func) ;
if (isOwned) partIntOwnedList->addOwned(*func) ;
compTermNorm.add(*norm,kFALSE) ;
}
}
// cout << "constructing special composite product" << endl ;
// cout << "compTermSet = " ; compTermSet.Print("1") ;
// WVE this doesn't work (yet) for some reason
//const char* name = makeRGPPName("SPECPROJ_",compTermSet,outerIntDeps,RooArgSet()) ;
//RooAbsReal* func = new RooGenProdProj(name,name,compTermSet,outerIntDeps,RooArgSet()) ;
//partIntList->add(*func) ;
//partIntOwnedList->addOwned(*func) ;
// WVE but this does, so we'll keep it for the moment
const char* prodname = makeRGPPName("SPECPROD_",compTermSet,outerIntDeps,RooArgSet(),isetRangeName) ;
RooProduct* prodtmp = new RooProduct(prodname,prodname,compTermSet) ;
partIntOwnedList->addOwned(*prodtmp) ;
const char* intname = makeRGPPName("SPECINT_",compTermSet,outerIntDeps,RooArgSet(),isetRangeName) ;
RooRealIntegral* inttmp = new RooRealIntegral(intname,intname,*prodtmp,outerIntDeps) ;
partIntOwnedList->addOwned(*inttmp) ;
partIntList->add(*inttmp);
partIntNormList->Add(compTermNorm.snapshot()) ;
delete tIter ;
}
}
delete gIter ;
// Store the partial integral list and return the assigned code ;
code = _partListMgr.setNormList(this,nset,iset,partIntList,RooNameReg::ptr(isetRangeName)) ;
_partOwnedListMgr.setNormList(this,nset,iset,partIntOwnedList,RooNameReg::ptr(isetRangeName)) ;
// Store the normalization set list in the cache using the index code of _partListMgr
_partNormListCache[code].Delete() ;
// We own content of 'norms' and transfer ownership to the cache
_partNormListCache[code] = *partIntNormList ;
// WVE PRINT SECTION
// cout << "RooProdPdf::getPIL(" << GetName() << "," << this << ") creating new configuration with code " << code << endl
// << " nset = " ; if (nset) nset->Print("1") ; else cout << "<none>" << endl ;
// cout << " iset" << (isetRangeName?(Form("[%s]",isetRangeName)):"") << " = " ; if (iset) iset->Print("1") ; else cout << "<none>" << endl ;
// cout << " Partial Integral List:" ;
// partIntList->Print("1") ;
// cout << " Partial Owned Integral List:" ;
// partIntOwnedList->Print("1") ;
// cout << " Normalization list: " << endl ; partIntNormList->Print("1") ;
// cout << endl ;
partList = partIntList ;
nsetList = _partNormListCache+code ;
// We own contents of all lists filled by factorizeProduct()
terms.Delete() ;
ints.Delete() ;
imp.Delete() ;
norms.Delete() ;
cross.Delete() ;
}
void RooProdPdf::groupProductTerms(RooLinkedList& groupedTerms, RooArgSet& outerIntDeps,
const RooLinkedList& terms, const RooLinkedList& norms,
const RooLinkedList& imps, const RooLinkedList& ints, const RooLinkedList& /*cross*/) const
{
// Start out with each term in its own group
TIterator* tIter = terms.MakeIterator() ;
RooArgSet* term ;
while((term=(RooArgSet*)tIter->Next())) {
RooLinkedList* group = new RooLinkedList ;
group->Add(term) ;
groupedTerms.Add(group) ;
}
delete tIter ;
// Make list of imported dependents that occur in any term
RooArgSet allImpDeps ;
TIterator* iIter = imps.MakeIterator() ;
RooArgSet *impDeps ;
while((impDeps=(RooArgSet*)iIter->Next())) {
allImpDeps.add(*impDeps,kFALSE) ;
}
delete iIter ;
// Make list of integrated dependents that occur in any term
RooArgSet allIntDeps ;
iIter = ints.MakeIterator() ;
RooArgSet *intDeps ;
while((intDeps=(RooArgSet*)iIter->Next())) {
allIntDeps.add(*intDeps,kFALSE) ;
}
delete iIter ;
// cout << "Complete lists of imported dependens" ; allImpDeps.Print("1") ;
// cout << "Complete lists of integrated dependens" ; allIntDeps.Print("1") ;
RooArgSet* tmp = (RooArgSet*) allIntDeps.selectCommon(allImpDeps) ;
outerIntDeps.removeAll() ;
outerIntDeps.add(*tmp) ;
delete tmp ;
// cout << "Integrated Dependents that need 'outer' integration: " ; outerIntDeps.Print("1") ;
// Now iteratively merge groups that should be (partially) integrated together
TIterator* oidIter = outerIntDeps.createIterator() ;
TIterator* glIter = groupedTerms.MakeIterator() ;
RooAbsArg* outerIntDep ;
while ((outerIntDep =(RooAbsArg*)oidIter->Next())) {
// cout << "merge for outerIntDep " << outerIntDep->GetName() << endl ;
// Collect groups that feature this dependent
RooLinkedList* newGroup = 0 ;
// Loop over groups
RooLinkedList* group ;
glIter->Reset() ;
Bool_t needMerge = kFALSE ;
while((group=(RooLinkedList*)glIter->Next())) {
// cout << "considering the following group:" << endl ;
// group->Print("1") ;
// cout << "----" << endl ;
// See if any term in this group depends in any ay on outerDepInt
RooArgSet* term ;
TIterator* tIter = group->MakeIterator() ;
while((term=(RooArgSet*)tIter->Next())) {
Int_t termIdx = terms.IndexOf(term) ;
RooArgSet* termNormDeps = (RooArgSet*) norms.At(termIdx) ;
RooArgSet* termIntDeps = (RooArgSet*) ints.At(termIdx) ;
RooArgSet* termImpDeps = (RooArgSet*) imps.At(termIdx) ;
if (termNormDeps->contains(*outerIntDep) ||
termIntDeps->contains(*outerIntDep) ||
termImpDeps->contains(*outerIntDep)) {
needMerge = kTRUE ;
}
}
if (needMerge) {
// cout << "group needs to be merged" << endl ;
// Create composite group if not yet existing
if (newGroup==0) {
newGroup = new RooLinkedList ;
}
// Add terms of this group to new term
tIter->Reset() ;
while((term=(RooArgSet*)tIter->Next())) {
// cout << "transferring group term to new merged group" << endl ;
newGroup->Add(term) ;
}
// Remove this group from list and delete it (but not its contents)
// cout << "removing group from groupList" << endl ;
groupedTerms.Remove(group) ;
delete group ;
}
delete tIter ;
}
// If a new group has been created to merge terms dependent on current outerIntDep, add it to group list
if (newGroup) {
// cout << "adding merged group to group list" << endl ;
groupedTerms.Add(newGroup) ;
// newGroup->Print("1") ;
// cout << "..." << endl ;
}
}
delete glIter ;
delete oidIter ;
}
RooAbsReal* RooProdPdf::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) const
{
// CASE I: factorizing term: term is integrated over all normalizing observables
// -----------------------------------------------------------------------------
// Check if all observbales of this term are integrated. If so the term cancels
if (termNSet.getSize()>0 && termNSet.getSize()==termISet.getSize() && isetRangeName==0) {
// Term factorizes
return 0 ;
}
// CASE II: Dropped terms: if term is entirely unnormalized, it should be dropped
// ------------------------------------------------------------------------------
if (nset && termNSet.getSize()==0) {
// Drop terms that are not asked to be normalized
return 0 ;
}
if (iset && termISet.getSize()>0) {
if (term->getSize()==1) {
// CASE IIIa: Normalized and partially integrated single PDF term
//---------------------------------------------------------------
TIterator* pIter = term->createIterator() ;
RooAbsPdf* pdf = (RooAbsPdf*) pIter->Next() ;
delete pIter ;
RooAbsReal* partInt = pdf->createIntegral(termISet,termNSet,isetRangeName) ;
partInt->setOperMode(operMode()) ;
isOwned=kTRUE ;
return partInt ;
} else {
// CASE IIIb: Normalized and partially integrated composite PDF term
//---------------------------------------------------------------
// Use auxiliary class RooGenProdProj to calculate this term
const char* name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
RooAbsReal* partInt = new RooGenProdProj(name,name,*term,termISet,termNSet,isetRangeName) ;
partInt->setOperMode(operMode()) ;
isOwned=kTRUE ;
return partInt ;
}
}
// CASE IVa: Normalized non-integrated composite PDF term
// -------------------------------------------------------
if (nset && nset->getSize()>0 && term->getSize()>1) {
// Composite term needs normalized integration
const char* name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
RooAbsReal* partInt = new RooGenProdProj(name,name,*term,termISet,termNSet,isetRangeName) ;
partInt->setOperMode(operMode()) ;
isOwned=kTRUE ;
return partInt ;
}
// CASE IVb: Normalized, non-integrated single PDF term
// -----------------------------------------------------
TIterator* pIter = term->createIterator() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)pIter->Next())) {
if (forceWrap) {
// Construct representative name of normalization wrapper
TString name(pdf->GetName()) ;
name.Append("_NORM[") ;
TIterator* nIter = termNSet.createIterator() ;
RooAbsArg* arg ;
Bool_t first(kTRUE) ;
while((arg=(RooAbsArg*)nIter->Next())) {
if (!first) {
name.Append(",") ;
} else {
first=kFALSE ;
}
name.Append(arg->GetName()) ;
}
name.Append("]") ;
delete nIter ;
RooAbsReal* partInt = new RooRealIntegral(name.Data(),name.Data(),*pdf,RooArgSet(),&termNSet) ;
isOwned=kTRUE ;
return partInt ;
} else {
isOwned=kFALSE ;
return pdf ;
}
}
delete pIter ;
cout << "RooProdPdf::processProductTerm(" << GetName() << ") unidentified term!!!" << endl ;
return 0;
}
const char* RooProdPdf::makeRGPPName(const char* pfx, const RooArgSet& term, const RooArgSet& iset, const RooArgSet& nset, const char* isetRangeName) const
{
// Make an appropriate name for a RooGenProdProj object in getPartIntList()
static TString pname ;
pname = pfx ;
TIterator* pIter = term.createIterator() ;
// Encode component names
Bool_t first(kTRUE) ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)pIter->Next())) {
if (first) {
first = kFALSE ;
} else {
pname.Append("_X_") ;
}
pname.Append(pdf->GetName()) ;
}
delete pIter ;
pname.Append(integralNameSuffix(iset,&nset,isetRangeName)) ;
return pname.Data() ;
}
Bool_t RooProdPdf::forceAnalyticalInt(const RooAbsArg& /*dep*/) const
{
return kTRUE ;
}
Int_t RooProdPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars,
const RooArgSet* normSet, const char* rangeName) const
{
// Determine which part (if any) of given integral can be performed analytically.
// If any analytical integration is possible, return integration scenario code.
//
// RooProdPdf implements two strategies in implementing analytical integrals
//
// First, PDF components whose entire set of dependents are requested to be integrated
// can be dropped from the product, as they will integrate out to 1 by construction
//:
// Second, RooProdPdf queries each remaining component PDF for its analytical integration
// capability of the requested set ('allVars'). It finds the largest common set of variables
// that can be integrated by all remaining components. If such a set exists, it reconfirms that
// each component is capable of analytically integrating the common set, and combines the components
// individual integration codes into a single integration code valid for RooProdPdf.
if (_forceNumInt) return 0 ;
// Declare that we can analytically integrate all requested observables
analVars.add(allVars) ;
// Retrieve (or create) the required partial integral list
Int_t code ;
RooArgList *plist ;
RooLinkedList *nlist ;
getPartIntList(normSet,&allVars,plist,nlist,code,rangeName) ;
return code+1 ;
}
Double_t RooProdPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
{
// Return analytical integral defined by given scenario code
// No integration scenario
if (code==0) {
return getVal(normSet) ;
}
// WVE needs adaptation for rangename feature
// Partial integration scenarios
RooArgList* partIntList = _partListMgr.getNormListByIndex(code-1) ;
// If cache has been sterilized, revive this slot
if (partIntList==0) {
RooArgSet* vars = getParameters(RooArgSet()) ;
RooArgSet* nset = _partListMgr.nameSet1ByIndex(code-1)->select(*vars) ;
RooArgSet* iset = _partListMgr.nameSet2ByIndex(code-1)->select(*vars) ;
Int_t code2(-1) ;
RooLinkedList *nlist ;
getPartIntList(nset,iset,partIntList,nlist,code2,rangeName) ;
delete vars ;
delete nset ;
delete iset ;
}
Double_t val = calculate(partIntList,_partNormListCache+code-1) ;
//cout << "RPP::aIWN(" << GetName() << ") value = " << val << endl ;
return val ;
}
Bool_t RooProdPdf::checkObservables(const RooArgSet* /*nset*/) const
{
// Check that none of the PDFs have overlapping dependents
return kFALSE ;
// Bool_t ret(kFALSE) ;
// _pdfIter->Reset() ;
// RooAbsPdf* pdf, *pdf2 ;
// TIterator* iter2 = _pdfList.createIterator() ;
// while((pdf = (RooAbsPdf*)_pdfIter->Next())) {
// *iter2 = *_pdfIter ;
// while((pdf2 = (RooAbsPdf*)iter2->Next())) {
// if (pdf->observableOverlaps(nset,*pdf2)) {
// cout << "RooProdPdf::checkObservables(" << GetName() << "): ERROR: PDFs " << pdf->GetName()
// << " and " << pdf2->GetName() << " have one or more dependents in common" << endl ;
// ret = kTRUE ;
// }
// }
// }
// delete iter2 ;
// return ret ;
}
RooAbsPdf::ExtendMode RooProdPdf::extendMode() const
{
return (_extendedIndex>=0) ? ((RooAbsPdf*)_pdfList.at(_extendedIndex))->extendMode() : CanNotBeExtended ;
}
Double_t RooProdPdf::expectedEvents(const RooArgSet* nset) const
{
assert(_extendedIndex>=0) ;
return ((RooAbsPdf*)_pdfList.at(_extendedIndex))->expectedEvents(nset) ;
}
RooAbsGenContext* RooProdPdf::genContext(const RooArgSet &vars, const RooDataSet *prototype,
const RooArgSet* auxProto, Bool_t verbose) const
{
if (_useDefaultGen) return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ;
return new RooProdGenContext(*this,vars,prototype,auxProto,verbose) ;
}
Int_t RooProdPdf::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
{
if (!_useDefaultGen) return 0 ;
// Find the subset directVars that only depend on a single PDF in the product
RooArgSet directSafe ;
TIterator* dIter = directVars.createIterator() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)dIter->Next())) {
if (isDirectGenSafe(*arg)) directSafe.add(*arg) ;
}
delete dIter ;
// Now find direct integrator for relevant components ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
Int_t code[64], n(0) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
RooArgSet pdfDirect ;
code[n] = pdf->getGenerator(directSafe,pdfDirect,staticInitOK) ;
if (code[n]!=0) {
generateVars.add(pdfDirect) ;
}
n++ ;
}
if (generateVars.getSize()>0) {
Int_t masterCode = _genCode.store(code,n) ;
return masterCode+1 ;
} else {
return 0 ;
}
}
void RooProdPdf::initGenerator(Int_t code)
{
if (!_useDefaultGen) return ;
const Int_t* codeList = _genCode.retrieve(code-1) ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
Int_t i(0) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
if (codeList[i]!=0) {
pdf->initGenerator(codeList[i]) ;
}
i++ ;
}
}
void RooProdPdf::generateEvent(Int_t code)
{
if (!_useDefaultGen) return ;
const Int_t* codeList = _genCode.retrieve(code-1) ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
Int_t i(0) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
if (codeList[i]!=0) {
pdf->generateEvent(codeList[i]) ;
}
i++ ;
}
}
void RooProdPdf::operModeHook()
{
Int_t i ;
for (i=0 ; i<_partListMgr.cacheSize() ; i++) {
RooArgList* plist = _partOwnedListMgr.getNormListByIndex(i) ;
if (plist) {
TIterator* iter = plist->createIterator() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)iter->Next())) {
arg->setOperMode(_operMode) ;
}
delete iter ;
}
}
return ;
}
void RooProdPdf::clearCache()
{
// cout << "RooProdPdf::clearCache(" << GetName() << ") sterilizing cache content" << endl ;
Int_t i ;
for (i=0 ; i<10 ; i++) {
_partNormListCache[i].Delete() ;
}
_partListMgr.sterilize() ;
_partOwnedListMgr.sterilize() ;
}
Bool_t RooProdPdf::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
{
// Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
clearCache() ;
return kFALSE ;
}
void RooProdPdf::printCompactTreeHook(ostream& os, const char* indent)
{
Int_t i ;
os << indent << "RooProdPdf begin partial integral cache" << endl ;
for (i=0 ; i<_partListMgr.cacheSize() ; i++) {
RooArgList* plist = _partListMgr.getNormListByIndex(i) ;
if (plist) {
TIterator* iter = plist->createIterator() ;
RooAbsArg* arg ;
TString indent2(indent) ;
indent2 += Form("[%d] ",i) ;
while((arg=(RooAbsArg*)iter->Next())) {
arg->printCompactTree(os,indent2) ;
}
delete iter ;
}
}
os << indent << "RooProdPdf end partial integral cache" << endl ;
}
Bool_t RooProdPdf::isDirectGenSafe(const RooAbsArg& arg) const
{
// Only override base class behaviour if default generator method is enabled
if (!_useDefaultGen) return RooAbsPdf::isDirectGenSafe(arg) ;
// Argument may appear in only one PDF component
_pdfIter->Reset() ;
RooAbsPdf* pdf, *thePdf(0) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
if (pdf->dependsOn(arg)) {
// Found PDF depending on arg
// If multiple PDFs depend on arg directGen is not safe
if (thePdf) return kFALSE ;
thePdf = pdf ;
}
}
// Forward call to relevant component PDF
return thePdf?(thePdf->isDirectGenSafe(arg)):kFALSE ;
}
RooArgSet* RooProdPdf::findPdfNSet(RooAbsPdf& pdf) const
{
// Look up user specified normalization set for given input PDF component
Int_t idx = _pdfList.index(&pdf) ;
if (idx<0) return 0 ;
return (RooArgSet*) _pdfNSetList.At(idx) ;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.