#include "RooFit.h"
#include "RooMsgService.h"
#include "TIterator.h"
#include "TIterator.h"
#include "TList.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooRealProxy.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooAddGenContext.h"
#include "RooRealConstant.h"
#include "RooNameReg.h"
#include "RooMsgService.h"
#include "RooRecursiveFraction.h"
#include "RooGlobalFunc.h"
#include "RooRealIntegral.h"
#include "RooTrace.h"
#include "Riostream.h"
#include <algorithm>
using namespace std;
ClassImp(RooAddPdf)
;
RooAddPdf::RooAddPdf() :
_refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
_refCoefRangeName(0),
_codeReg(10),
_snormList(0),
_recursive(kFALSE)
{
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_coefCache = new Double_t[100] ;
_coefErrCount = _errorCount ;
TRACE_CREATE
}
RooAddPdf::RooAddPdf(const char *name, const char *title) :
RooAbsPdf(name,title),
_refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
_refCoefRangeName(0),
_projectCoefs(kFALSE),
_projCacheMgr(this,10),
_codeReg(10),
_pdfList("!pdfs","List of PDFs",this),
_coefList("!coefficients","List of coefficients",this),
_snormList(0),
_haveLastCoef(kFALSE),
_allExtendable(kFALSE),
_recursive(kFALSE)
{
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_coefCache = new Double_t[100] ;
_coefErrCount = _errorCount ;
TRACE_CREATE
}
RooAddPdf::RooAddPdf(const char *name, const char *title,
RooAbsPdf& pdf1, RooAbsPdf& pdf2, RooAbsReal& coef1) :
RooAbsPdf(name,title),
_refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
_refCoefRangeName(0),
_projectCoefs(kFALSE),
_projCacheMgr(this,10),
_codeReg(10),
_pdfList("!pdfs","List of PDFs",this),
_coefList("!coefficients","List of coefficients",this),
_haveLastCoef(kFALSE),
_allExtendable(kFALSE),
_recursive(kFALSE)
{
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_pdfList.add(pdf1) ;
_pdfList.add(pdf2) ;
_coefList.add(coef1) ;
_coefCache = new Double_t[_pdfList.getSize()] ;
_coefErrCount = _errorCount ;
TRACE_CREATE
}
RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList, const RooArgList& inCoefList, Bool_t recursiveFractions) :
RooAbsPdf(name,title),
_refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
_refCoefRangeName(0),
_projectCoefs(kFALSE),
_projCacheMgr(this,10),
_codeReg(10),
_pdfList("!pdfs","List of PDFs",this),
_coefList("!coefficients","List of coefficients",this),
_haveLastCoef(kFALSE),
_allExtendable(kFALSE),
_recursive(kFALSE)
{
if (inPdfList.getSize()>inCoefList.getSize()+1 || inPdfList.getSize()<inCoefList.getSize()) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName()
<< ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
assert(0) ;
}
if (recursiveFractions && inPdfList.getSize()!=inCoefList.getSize()+1) {
coutW(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName()
<< ") WARNING inconsistent input: recursive fractions options can only be used if Npdf=Ncoef+1, ignoring recursive fraction setting" << endl ;
}
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
TIterator* pdfIter = inPdfList.createIterator() ;
TIterator* coefIter = inCoefList.createIterator() ;
RooAbsPdf* pdf ;
RooAbsReal* coef ;
RooArgList partinCoefList ;
Bool_t first(kTRUE) ;
while((coef = (RooAbsPdf*)coefIter->Next())) {
pdf = (RooAbsPdf*) pdfIter->Next() ;
if (!pdf) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName()
<< ") number of pdfs and coefficients inconsistent, must have Npdf=Ncoef or Npdf=Ncoef+1" << endl ;
assert(0) ;
}
if (!dynamic_cast<RooAbsReal*>(coef)) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal, ignored" << endl ;
continue ;
}
if (!dynamic_cast<RooAbsReal*>(pdf)) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName() << " is not of type RooAbsPdf, ignored" << endl ;
continue ;
}
_pdfList.add(*pdf) ;
if (recursiveFractions) {
partinCoefList.add(*coef) ;
if (first) {
first = kFALSE ;
_coefList.add(*coef) ;
} else {
RooAbsReal* rfrac = new RooRecursiveFraction(Form("%s_recursive_fraction_%s",GetName(),pdf->GetName()),"Recursive Fraction",partinCoefList) ;
addOwnedComponents(*rfrac) ;
_coefList.add(*rfrac) ;
}
} else {
_coefList.add(*coef) ;
}
}
pdf = (RooAbsPdf*) pdfIter->Next() ;
if (pdf) {
if (!dynamic_cast<RooAbsReal*>(pdf)) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") last pdf " << coef->GetName() << " is not of type RooAbsPdf, fatal error" << endl ;
assert(0) ;
}
_pdfList.add(*pdf) ;
if (recursiveFractions) {
partinCoefList.add(RooFit::RooConst(1)) ;
RooAbsReal* rfrac = new RooRecursiveFraction(Form("%s_recursive_fraction_%s",GetName(),pdf->GetName()),"Recursive Fraction",partinCoefList) ;
addOwnedComponents(*rfrac) ;
_coefList.add(*rfrac) ;
_haveLastCoef=kTRUE ;
}
} else {
_haveLastCoef=kTRUE ;
}
delete pdfIter ;
delete coefIter ;
_coefCache = new Double_t[_pdfList.getSize()] ;
_coefErrCount = _errorCount ;
_recursive = recursiveFractions ;
TRACE_CREATE
}
RooAddPdf::RooAddPdf(const char *name, const char *title, const RooArgList& inPdfList) :
RooAbsPdf(name,title),
_refCoefNorm("!refCoefNorm","Reference coefficient normalization set",this,kFALSE,kFALSE),
_refCoefRangeName(0),
_projectCoefs(kFALSE),
_projCacheMgr(this,10),
_pdfList("!pdfs","List of PDFs",this),
_coefList("!coefficients","List of coefficients",this),
_haveLastCoef(kFALSE),
_allExtendable(kTRUE),
_recursive(kFALSE)
{
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
TIterator* pdfIter = inPdfList.createIterator() ;
RooAbsPdf* pdf ;
while((pdf = (RooAbsPdf*) pdfIter->Next())) {
if (!dynamic_cast<RooAbsReal*>(pdf)) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName() << " is not of type RooAbsPdf, ignored" << endl ;
continue ;
}
if (!pdf->canBeExtended()) {
coutE(InputArguments) << "RooAddPdf::RooAddPdf(" << GetName() << ") pdf " << pdf->GetName() << " is not extendable, ignored" << endl ;
continue ;
}
_pdfList.add(*pdf) ;
}
delete pdfIter ;
_coefCache = new Double_t[_pdfList.getSize()] ;
_coefErrCount = _errorCount ;
TRACE_CREATE
}
RooAddPdf::RooAddPdf(const RooAddPdf& other, const char* name) :
RooAbsPdf(other,name),
_refCoefNorm("!refCoefNorm",this,other._refCoefNorm),
_refCoefRangeName((TNamed*)other._refCoefRangeName),
_projectCoefs(other._projectCoefs),
_projCacheMgr(other._projCacheMgr,this),
_codeReg(other._codeReg),
_pdfList("!pdfs",this,other._pdfList),
_coefList("!coefficients",this,other._coefList),
_haveLastCoef(other._haveLastCoef),
_allExtendable(other._allExtendable),
_recursive(other._recursive)
{
_pdfIter = _pdfList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_coefCache = new Double_t[_pdfList.getSize()] ;
_coefErrCount = _errorCount ;
TRACE_CREATE
}
RooAddPdf::~RooAddPdf()
{
delete _pdfIter ;
delete _coefIter ;
if (_coefCache) delete[] _coefCache ;
TRACE_DESTROY
}
void RooAddPdf::fixCoefNormalization(const RooArgSet& refCoefNorm)
{
if (refCoefNorm.getSize()==0) {
_projectCoefs = kFALSE ;
return ;
}
_projectCoefs = kTRUE ;
_refCoefNorm.removeAll() ;
_refCoefNorm.add(refCoefNorm) ;
_projCacheMgr.reset() ;
}
void RooAddPdf::fixCoefRange(const char* rangeName)
{
_refCoefRangeName = (TNamed*)RooNameReg::ptr(rangeName) ;
if (_refCoefRangeName) _projectCoefs = kTRUE ;
}
RooAddPdf::CacheElem* RooAddPdf::getProjCache(const RooArgSet* nset, const RooArgSet* iset, const char* rangeName) const
{
CacheElem* cache = (CacheElem*) _projCacheMgr.getObj(nset,iset,0,rangeName) ;
if (cache) {
return cache ;
}
cache = new CacheElem ;
RooArgSet *fullDepList = getObservables(nset) ;
if (iset) {
fullDepList->remove(*iset,kTRUE,kTRUE) ;
}
_pdfIter->Reset() ;
_coefIter->Reset() ;
RooAbsPdf* pdf ;
RooAbsReal* coef ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
coef=(RooAbsPdf*)_coefIter->Next() ;
RooArgSet supNSet(*fullDepList) ;
RooArgSet* pdfDeps = pdf->getObservables(nset) ;
if (pdfDeps) {
supNSet.remove(*pdfDeps,kTRUE,kTRUE) ;
delete pdfDeps ;
}
RooArgSet* coefDeps = coef ? coef->getObservables(nset) : 0 ;
if (coefDeps) {
supNSet.remove(*coefDeps,kTRUE,kTRUE) ;
delete coefDeps ;
}
RooAbsReal* snorm ;
TString name(GetName()) ;
name.Append("_") ;
name.Append(pdf->GetName()) ;
name.Append("_SupNorm") ;
cache->_needSupNorm = kFALSE ;
if (supNSet.getSize()>0) {
snorm = new RooRealIntegral(name,"Supplemental normalization integral",RooRealConstant::value(1.0),supNSet) ;
cxcoutD(Caching) << "RooAddPdf " << GetName() << " making supplemental normalization set " << supNSet << " for pdf component " << pdf->GetName() << endl ;
cache->_needSupNorm = kTRUE ;
} else {
snorm = new RooRealVar(name,"Unit Supplemental normalization integral",1.0) ;
}
cache->_suppNormList.addOwned(*snorm) ;
}
delete fullDepList ;
if (_verboseEval>1) {
cxcoutD(Caching) << "RooAddPdf::syncSuppNormList(" << GetName() << ") synching supplemental normalization list for norm" << (nset?*nset:RooArgSet()) << endl ;
if dologD(Caching) {
cache->_suppNormList.Print("v") ;
}
}
if (!_projectCoefs && !rangeName) {
_projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
return cache ;
}
RooArgSet* nset2 = nset ? getObservables(nset) : new RooArgSet() ;
cxcoutD(Caching) << "RooAddPdf(" << GetName() << ")::getPC nset = " << (nset?*nset:RooArgSet()) << " nset2 = " << *nset2 << endl ;
if (nset2->getSize()==0 && _refCoefNorm.getSize()!=0) {
nset2->add(_refCoefNorm) ;
if (_refCoefRangeName) {
rangeName = RooNameReg::str(_refCoefRangeName) ;
}
}
if (!nset2->equals(_refCoefNorm) || _refCoefRangeName !=0 || rangeName !=0 || _normRange.Length()>0) {
cxcoutD(Caching) << "ALEX: RooAddPdf::syncCoefProjList(" << GetName() << ") projecting coefficients from "
<< *nset2 << (rangeName?":":"") << (rangeName?rangeName:"")
<< " to " << ((_refCoefNorm.getSize()>0)?_refCoefNorm:*nset2) << (_refCoefRangeName?":":"") << (_refCoefRangeName?RooNameReg::str(_refCoefRangeName):"") << endl ;
_pdfIter->Reset() ;
RooAbsPdf* thePdf ;
while((thePdf=(RooAbsPdf*)_pdfIter->Next())) {
RooAbsReal* pdfProj ;
if (!nset2->equals(_refCoefNorm)) {
pdfProj = thePdf->createIntegral(*nset2,_refCoefNorm,_normRange.Length()>0?_normRange.Data():0) ;
pdfProj->setOperMode(operMode()) ;
cxcoutD(Caching) << "RooAddPdf(" << GetName() << ")::getPC nset2(" << *nset2 << ")!=_refCoefNorm(" << _refCoefNorm << ") --> pdfProj = " << pdfProj->GetName() << endl ;
} else {
TString name(GetName()) ;
name.Append("_") ;
name.Append(thePdf->GetName()) ;
name.Append("_ProjectNorm") ;
pdfProj = new RooRealVar(name,"Unit Projection normalization integral",1.0) ;
cxcoutD(Caching) << "RooAddPdf(" << GetName() << ")::getPC nset2(" << *nset2 << ")==_refCoefNorm(" << _refCoefNorm << ") --> pdfProj = " << pdfProj->GetName() << endl ;
}
cache->_projList.addOwned(*pdfProj) ;
cxcoutD(Caching) << " RooAddPdf::syncCoefProjList(" << GetName() << ") PP = " << pdfProj->GetName() << endl ;
RooArgSet supNormSet(_refCoefNorm) ;
RooArgSet* deps = thePdf->getParameters(RooArgSet()) ;
supNormSet.remove(*deps,kTRUE,kTRUE) ;
delete deps ;
RooAbsReal* snorm ;
TString name(GetName()) ;
name.Append("_") ;
name.Append(thePdf->GetName()) ;
name.Append("_ProjSupNorm") ;
if (supNormSet.getSize()>0 && !nset2->equals(_refCoefNorm) ) {
snorm = new RooRealIntegral(name,"Projection Supplemental normalization integral",
RooRealConstant::value(1.0),supNormSet) ;
} else {
snorm = new RooRealVar(name,"Unit Projection Supplemental normalization integral",1.0) ;
}
cxcoutD(Caching) << " RooAddPdf::syncCoefProjList(" << GetName() << ") SN = " << snorm->GetName() << endl ;
cache->_suppProjList.addOwned(*snorm) ;
RooAbsReal* rangeProj1 ;
RooArgSet* tmpObs = thePdf->getObservables(_refCoefNorm) ;
RooAbsArg* obsArg ;
TIterator* iter = tmpObs->createIterator() ;
Bool_t allIdent = kTRUE ;
while((obsArg=(RooAbsArg*)iter->Next())) {
RooRealVar* rvarg = dynamic_cast<RooRealVar*>(obsArg) ;
if (rvarg) {
if (rvarg->getMin(RooNameReg::str(_refCoefRangeName))!=rvarg->getMin() ||
rvarg->getMax(RooNameReg::str(_refCoefRangeName))!=rvarg->getMax()) {
allIdent=kFALSE ;
}
}
}
delete iter ;
delete tmpObs ;
if (_refCoefRangeName && _refCoefNorm.getSize()>0 && !allIdent) {
RooArgSet* tmp = thePdf->getObservables(_refCoefNorm) ;
rangeProj1 = thePdf->createIntegral(*tmp,*tmp,RooNameReg::str(_refCoefRangeName)) ;
delete tmp ;
} else {
TString theName(GetName()) ;
theName.Append("_") ;
theName.Append(thePdf->GetName()) ;
theName.Append("_RangeNorm1") ;
rangeProj1 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
}
cxcoutD(Caching) << " RooAddPdf::syncCoefProjList(" << GetName() << ") R1 = " << rangeProj1->GetName() << endl ;
cache->_refRangeProjList.addOwned(*rangeProj1) ;
RooAbsReal* rangeProj2 ;
cxcoutD(Caching) << "RooAddPdf::syncCoefProjList(" << GetName() << ") rangename = " << (rangeName?rangeName:"<null>")
<< " nset = " << (nset?*nset:RooArgSet()) << endl ;
if (rangeName && _refCoefNorm.getSize()>0) {
rangeProj2 = thePdf->createIntegral(_refCoefNorm,_refCoefNorm,rangeName) ;
} else if (_normRange.Length()>0) {
RooArgSet* tmp = thePdf->getObservables(_refCoefNorm) ;
rangeProj2 = thePdf->createIntegral(*tmp,*tmp,_normRange.Data()) ;
delete tmp ;
} else {
TString theName(GetName()) ;
theName.Append("_") ;
theName.Append(thePdf->GetName()) ;
theName.Append("_RangeNorm2") ;
rangeProj2 = new RooRealVar(theName,"Unit range normalization integral",1.0) ;
}
cxcoutD(Caching) << " RooAddPdf::syncCoefProjList(" << GetName() << ") R2 = " << rangeProj2->GetName() << endl ;
cache->_rangeProjList.addOwned(*rangeProj2) ;
}
}
delete nset2 ;
_projCacheMgr.setObj(nset,iset,cache,RooNameReg::ptr(rangeName)) ;
return cache ;
}
void RooAddPdf::updateCoefficients(CacheElem& cache, const RooArgSet* nset) const
{
Int_t i ;
if (_allExtendable) {
Double_t coefSum(0) ;
RooFIter it=_pdfList.fwdIterator() ; i=0 ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)it.next())) {
_coefCache[i] = pdf->expectedEvents(_refCoefNorm.getSize()>0?&_refCoefNorm:nset) ;
coefSum += _coefCache[i] ;
i++ ;
}
if (coefSum==0.) {
coutW(Eval) << "RooAddPdf::updateCoefCache(" << GetName() << ") WARNING: total number of expected events is 0" << endl ;
} else {
Int_t siz = _pdfList.getSize() ;
for (i=0 ; i<siz ; i++) {
_coefCache[i] /= coefSum ;
}
}
} else {
if (_haveLastCoef) {
Double_t coefSum(0) ;
RooFIter it=_coefList.fwdIterator() ; i=0 ;
RooAbsReal* coef ;
while((coef=(RooAbsReal*)it.next())) {
_coefCache[i] = coef->getVal(nset) ;
coefSum += _coefCache[i] ;
i++ ;
}
if (coefSum==0.) {
coutW(Eval) << "RooAddPdf::updateCoefCache(" << GetName() << ") WARNING: sum of coefficients is zero 0" << endl ;
} else {
Int_t siz = _coefList.getSize() ;
for (i=0 ; i<siz ; i++) {
_coefCache[i] /= coefSum ;
}
}
} else {
Double_t lastCoef(1) ;
RooFIter it=_coefList.fwdIterator() ; i=0 ;
RooAbsReal* coef ;
while((coef=(RooAbsReal*)it.next())) {
_coefCache[i] = coef->getVal(nset) ;
lastCoef -= _coefCache[i] ;
i++ ;
}
_coefCache[_coefList.getSize()] = lastCoef ;
if ((lastCoef<-1e-05 || (lastCoef-1)>1e-5) && _coefErrCount-->0) {
coutW(Eval) << "RooAddPdf::updateCoefCache(" << GetName()
<< " WARNING: sum of PDF coefficients not in range [0-1], value="
<< 1-lastCoef ;
if (_coefErrCount==0) {
coutW(Eval) << " (no more will be printed)" ;
}
coutW(Eval) << endl ;
}
}
}
if ((!_projectCoefs && _normRange.Length()==0) || cache._projList.getSize()==0) {
return ;
}
Double_t coefSum(0) ;
for (i=0 ; i<_pdfList.getSize() ; i++) {
Bool_t _tmp = _globalSelectComp ;
RooAbsPdf::globalSelectComp(kTRUE) ;
RooAbsReal* pp = ((RooAbsReal*)cache._projList.at(i)) ;
RooAbsReal* sn = ((RooAbsReal*)cache._suppProjList.at(i)) ;
RooAbsReal* r1 = ((RooAbsReal*)cache._refRangeProjList.at(i)) ;
RooAbsReal* r2 = ((RooAbsReal*)cache._rangeProjList.at(i)) ;
Double_t proj = pp->getVal()/sn->getVal()*(r2->getVal()/r1->getVal()) ;
RooAbsPdf::globalSelectComp(_tmp) ;
_coefCache[i] *= proj ;
coefSum += _coefCache[i] ;
}
for (i=0 ; i<_pdfList.getSize() ; i++) {
_coefCache[i] /= coefSum ;
cxcoutD(Caching) << " ALEX: POST-SYNC coef[" << i << "] = " << _coefCache[i]
<< " ( _coefCache[i]/coefSum = " << _coefCache[i]*coefSum << "/" << coefSum << " ) "<< endl ;
}
}
Double_t RooAddPdf::evaluate() const
{
const RooArgSet* nset = _normSet ;
if (nset==0 || nset->getSize()==0) {
if (_refCoefNorm.getSize()!=0) {
nset = &_refCoefNorm ;
}
}
CacheElem* cache = getProjCache(nset) ;
updateCoefficients(*cache,nset) ;
RooAbsPdf* pdf ;
Double_t value(0) ;
Int_t i(0) ;
RooFIter pi = _pdfList.fwdIterator() ;
if (cache->_needSupNorm) {
Double_t snormVal ;
while((pdf = (RooAbsPdf*)pi.next())) {
snormVal = ((RooAbsReal*)cache->_suppNormList.at(i))->getVal() ;
Double_t pdfVal = pdf->getVal(nset) ;
if (pdf->isSelectedComp()) {
value += pdfVal*_coefCache[i]/snormVal ;
}
i++ ;
}
} else {
while((pdf = (RooAbsPdf*)pi.next())) {
Double_t pdfVal = pdf->getVal(nset) ;
if (pdf->isSelectedComp()) {
value += pdfVal*_coefCache[i] ;
}
i++ ;
}
}
return value ;
}
void RooAddPdf::resetErrorCounters(Int_t resetValue)
{
RooAbsPdf::resetErrorCounters(resetValue) ;
_coefErrCount = resetValue ;
}
Bool_t RooAddPdf::checkObservables(const RooArgSet* nset) const
{
Bool_t ret(kFALSE) ;
_pdfIter->Reset() ;
_coefIter->Reset() ;
RooAbsReal* coef ;
RooAbsReal* pdf ;
while((coef=(RooAbsReal*)_coefIter->Next())) {
pdf = (RooAbsReal*)_pdfIter->Next() ;
if (pdf->observableOverlaps(nset,*coef)) {
coutE(InputArguments) << "RooAddPdf::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName()
<< " and PDF " << pdf->GetName() << " have one or more dependents in common" << endl ;
ret = kTRUE ;
}
}
return ret ;
}
Int_t RooAddPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars,
const RooArgSet* normSet, const char* rangeName) const
{
RooArgSet* allDepVars = getObservables(allVars) ;
RooArgSet allAnalVars(*allDepVars) ;
delete allDepVars ;
TIterator* avIter = allVars.createIterator() ;
Int_t n(0) ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
RooArgSet subAnalVars ;
pdf->getAnalyticalIntegralWN(allVars,subAnalVars,normSet,rangeName) ;
avIter->Reset() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)avIter->Next())) {
if (!subAnalVars.find(arg->GetName()) && pdf->dependsOn(*arg)) {
allAnalVars.remove(*arg,kTRUE,kTRUE) ;
}
}
n++ ;
}
if (allAnalVars.getSize()==0) {
delete avIter ;
return 0 ;
}
_pdfIter->Reset() ;
n=0 ;
std::vector<Int_t> subCode(_pdfList.getSize());
Bool_t allOK(kTRUE) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
RooArgSet subAnalVars ;
RooArgSet* allAnalVars2 = pdf->getObservables(allAnalVars) ;
subCode[n] = pdf->getAnalyticalIntegralWN(*allAnalVars2,subAnalVars,normSet,rangeName) ;
if (subCode[n]==0 && allAnalVars2->getSize()>0) {
coutE(InputArguments) << "RooAddPdf::getAnalyticalIntegral(" << GetName() << ") WARNING: component PDF " << pdf->GetName()
<< " advertises inconsistent set of integrals (e.g. (X,Y) but not X or Y individually."
<< " Distributed analytical integration disabled. Please fix PDF" << endl ;
allOK = kFALSE ;
}
delete allAnalVars2 ;
n++ ;
}
if (!allOK) {
delete avIter ;
return 0 ;
}
analVars.add(allAnalVars) ;
RooArgSet* intSet = new RooArgSet(allAnalVars) ;
Int_t masterCode = _codeReg.store(subCode,intSet)+1 ;
delete avIter ;
return masterCode ;
}
Double_t RooAddPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
{
if (code==0) {
return getVal(normSet) ;
}
RooArgSet* intSet ;
const std::vector<Int_t>& subCode = _codeReg.retrieve(code-1,intSet) ;
if (subCode.empty()) {
coutE(InputArguments) << "RooAddPdf::analyticalIntegral(" << GetName() << "): ERROR unrecognized integration code, " << code << endl ;
assert(0) ;
}
cxcoutD(Caching) << "RooAddPdf::aiWN(" << GetName() << ") calling getProjCache with nset = " << (normSet?*normSet:RooArgSet()) << endl ;
if ((normSet==0 || normSet->getSize()==0) && _refCoefNorm.getSize()>0) {
normSet = &_refCoefNorm ;
}
CacheElem* cache = getProjCache(normSet,intSet,0) ;
updateCoefficients(*cache,normSet) ;
Double_t value(0) ;
_pdfIter->Reset() ;
_coefIter->Reset() ;
RooAbsPdf* pdf ;
Double_t snormVal ;
Int_t i(0) ;
RooArgList* snormSet = (cache->_suppNormList.getSize()>0) ? &cache->_suppNormList : 0 ;
while((pdf = (RooAbsPdf*)_pdfIter->Next())) {
if (_coefCache[i]) {
snormVal = snormSet ? ((RooAbsReal*) cache->_suppNormList.at(i))->getVal() : 1.0 ;
Double_t val = pdf->analyticalIntegralWN(subCode[i],normSet,rangeName) ;
if (pdf->isSelectedComp()) {
value += val*_coefCache[i]/snormVal ;
}
}
i++ ;
}
return value ;
}
Double_t RooAddPdf::expectedEvents(const RooArgSet* nset) const
{
Double_t expectedTotal(0.0);
cxcoutD(Caching) << "RooAddPdf::expectedEvents(" << GetName() << ") calling getProjCache with nset = " << (nset?*nset:RooArgSet()) << endl ;
CacheElem* cache = getProjCache(nset) ;
updateCoefficients(*cache,nset) ;
if (cache->_rangeProjList.getSize()>0) {
RooFIter iter1 = cache->_refRangeProjList.fwdIterator() ;
RooFIter iter2 = cache->_rangeProjList.fwdIterator() ;
RooFIter iter3 = _pdfList.fwdIterator() ;
if (_allExtendable) {
RooAbsPdf* pdf ;
while ((pdf=(RooAbsPdf*)iter3.next())) {
RooAbsReal* r1 = (RooAbsReal*)iter1.next() ;
RooAbsReal* r2 = (RooAbsReal*)iter2.next() ;
expectedTotal += (r2->getVal()/r1->getVal()) * pdf->expectedEvents(nset) ;
}
} else {
RooFIter citer = _coefList.fwdIterator() ;
RooAbsReal* coef ;
while((coef=(RooAbsReal*)citer.next())) {
Double_t ncomp = coef->getVal(nset) ;
RooAbsReal* r1 = (RooAbsReal*)iter1.next() ;
RooAbsReal* r2 = (RooAbsReal*)iter2.next() ;
expectedTotal += (r2->getVal()/r1->getVal()) * ncomp ;
}
}
} else {
if (_allExtendable) {
RooFIter iter = _pdfList.fwdIterator() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)iter.next())) {
expectedTotal += pdf->expectedEvents(nset) ;
}
} else {
RooFIter citer = _coefList.fwdIterator() ;
RooAbsReal* coef ;
while((coef=(RooAbsReal*)citer.next())) {
Double_t ncomp = coef->getVal(nset) ;
expectedTotal += ncomp ;
}
}
}
return expectedTotal ;
}
void RooAddPdf::selectNormalization(const RooArgSet* depSet, Bool_t force)
{
if (!force && _refCoefNorm.getSize()!=0) {
return ;
}
if (!depSet) {
fixCoefNormalization(RooArgSet()) ;
return ;
}
RooArgSet* myDepSet = getObservables(depSet) ;
fixCoefNormalization(*myDepSet) ;
delete myDepSet ;
}
void RooAddPdf::selectNormalizationRange(const char* rangeName, Bool_t force)
{
if (!force && _refCoefRangeName) {
return ;
}
fixCoefRange(rangeName) ;
}
RooAbsGenContext* RooAddPdf::genContext(const RooArgSet &vars, const RooDataSet *prototype,
const RooArgSet* auxProto, Bool_t verbose) const
{
return new RooAddGenContext(*this,vars,prototype,auxProto,verbose) ;
}
RooArgList RooAddPdf::CacheElem::containedArgs(Action)
{
RooArgList allNodes;
allNodes.add(_projList) ;
allNodes.add(_suppProjList) ;
allNodes.add(_refRangeProjList) ;
allNodes.add(_rangeProjList) ;
return allNodes ;
}
std::list<Double_t>* RooAddPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
list<Double_t>* sumHint = 0 ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
Bool_t needClean(kFALSE) ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
list<Double_t>* pdfHint = pdf->plotSamplingHint(obs,xlo,xhi) ;
if (pdfHint) {
if (!sumHint) {
sumHint = pdfHint ;
} else {
list<Double_t>* newSumHint = new list<Double_t>(sumHint->size()+pdfHint->size()) ;
merge(pdfHint->begin(),pdfHint->end(),sumHint->begin(),sumHint->end(),newSumHint->begin()) ;
delete sumHint ;
sumHint = newSumHint ;
needClean = kTRUE ;
}
}
}
if (needClean) {
list<Double_t>::iterator new_end = unique(sumHint->begin(),sumHint->end()) ;
sumHint->erase(new_end,sumHint->end()) ;
}
return sumHint ;
}
std::list<Double_t>* RooAddPdf::binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
list<Double_t>* sumBinB = 0 ;
Bool_t needClean(kFALSE) ;
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
list<Double_t>* pdfBinB = pdf->binBoundaries(obs,xlo,xhi) ;
if (pdfBinB) {
if (!sumBinB) {
sumBinB = pdfBinB ;
} else {
list<Double_t>* newSumBinB = new list<Double_t>(sumBinB->size()+pdfBinB->size()) ;
merge(pdfBinB->begin(),pdfBinB->end(),sumBinB->begin(),sumBinB->end(),newSumBinB->begin()) ;
delete sumBinB ;
delete pdfBinB ;
sumBinB = newSumBinB ;
needClean = kTRUE ;
}
}
}
if (needClean) {
list<Double_t>::iterator new_end = unique(sumBinB->begin(),sumBinB->end()) ;
sumBinB->erase(new_end,sumBinB->end()) ;
}
return sumBinB ;
}
Bool_t RooAddPdf::isBinnedDistribution(const RooArgSet& obs) const
{
_pdfIter->Reset() ;
RooAbsPdf* pdf ;
while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
if (pdf->dependsOn(obs) && !pdf->isBinnedDistribution(obs)) {
return kFALSE ;
}
}
return kTRUE ;
}
void RooAddPdf::setCacheAndTrackHints(RooArgSet& trackNodes)
{
RooFIter aiter = pdfList().fwdIterator() ;
RooAbsArg* aarg ;
while ((aarg=aiter.next())) {
if (aarg->canNodeBeCached()==Always) {
trackNodes.add(*aarg) ;
}
}
}
void RooAddPdf::printMetaArgs(ostream& os) const
{
_pdfIter->Reset() ;
_coefIter->Reset() ;
Bool_t first(kTRUE) ;
RooAbsArg* coef, *pdf ;
if (_coefList.getSize()!=0) {
while((coef=(RooAbsArg*)_coefIter->Next())) {
if (!first) {
os << " + " ;
} else {
first = kFALSE ;
}
pdf=(RooAbsArg*)_pdfIter->Next() ;
os << coef->GetName() << " * " << pdf->GetName() ;
}
pdf = (RooAbsArg*) _pdfIter->Next() ;
if (pdf) {
os << " + [%] * " << pdf->GetName() ;
}
} else {
while((pdf=(RooAbsArg*)_pdfIter->Next())) {
if (!first) {
os << " + " ;
} else {
first = kFALSE ;
}
os << pdf->GetName() ;
}
}
os << " " ;
}