#include "RooFit.h"
#include "Riostream.h"
#include "TIterator.h"
#include "TList.h"
#include "RooRealSumPdf.h"
#include "RooRealProxy.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooAddGenContext.h"
#include "RooRealConstant.h"
#include "RooRealIntegral.h"
#include "RooMsgService.h"
#include "RooNameReg.h"
#include <memory>
#include <algorithm>
#include "TError.h"
using namespace std;
ClassImp(RooRealSumPdf)
;
Bool_t RooRealSumPdf::_doFloorGlobal = kFALSE ;
RooRealSumPdf::RooRealSumPdf()
{
_funcIter = _funcList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_extended = kFALSE ;
_doFloor = kFALSE ;
}
RooRealSumPdf::RooRealSumPdf(const char *name, const char *title) :
RooAbsPdf(name,title),
_normIntMgr(this,10),
_haveLastCoef(kFALSE),
_funcList("!funcList","List of functions",this),
_coefList("!coefList","List of coefficients",this),
_extended(kFALSE),
_doFloor(kFALSE)
{
_funcIter = _funcList.createIterator() ;
_coefIter = _coefList.createIterator() ;
}
RooRealSumPdf::RooRealSumPdf(const char *name, const char *title,
RooAbsReal& func1, RooAbsReal& func2, RooAbsReal& coef1) :
RooAbsPdf(name,title),
_normIntMgr(this,10),
_haveLastCoef(kFALSE),
_funcList("!funcList","List of functions",this),
_coefList("!coefList","List of coefficients",this),
_extended(kFALSE),
_doFloor(kFALSE)
{
_funcIter = _funcList.createIterator() ;
_coefIter = _coefList.createIterator() ;
_funcList.add(func1) ;
_funcList.add(func2) ;
_coefList.add(coef1) ;
}
RooRealSumPdf::RooRealSumPdf(const char *name, const char *title, const RooArgList& inFuncList, const RooArgList& inCoefList, Bool_t extended) :
RooAbsPdf(name,title),
_normIntMgr(this,10),
_haveLastCoef(kFALSE),
_funcList("!funcList","List of functions",this),
_coefList("!coefList","List of coefficients",this),
_extended(extended),
_doFloor(kFALSE)
{
if (!(inFuncList.getSize()==inCoefList.getSize()+1 || inFuncList.getSize()==inCoefList.getSize())) {
coutE(InputArguments) << "RooRealSumPdf::RooRealSumPdf(" << GetName()
<< ") number of pdfs and coefficients inconsistent, must have Nfunc=Ncoef or Nfunc=Ncoef+1" << endl ;
assert(0) ;
}
_funcIter = _funcList.createIterator() ;
_coefIter = _coefList.createIterator() ;
TIterator* funcIter = inFuncList.createIterator() ;
TIterator* coefIter = inCoefList.createIterator() ;
RooAbsArg* func ;
RooAbsArg* coef ;
while((coef = (RooAbsArg*)coefIter->Next())) {
func = (RooAbsArg*) funcIter->Next() ;
if (!dynamic_cast<RooAbsReal*>(coef)) {
coutW(InputArguments) << "RooRealSumPdf::RooRealSumPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal, ignored" << endl ;
continue ;
}
if (!dynamic_cast<RooAbsReal*>(func)) {
coutW(InputArguments) << "RooRealSumPdf::RooRealSumPdf(" << GetName() << ") func " << func->GetName() << " is not of type RooAbsReal, ignored" << endl ;
continue ;
}
_funcList.add(*func) ;
_coefList.add(*coef) ;
}
func = (RooAbsReal*) funcIter->Next() ;
if (func) {
if (!dynamic_cast<RooAbsReal*>(func)) {
coutE(InputArguments) << "RooRealSumPdf::RooRealSumPdf(" << GetName() << ") last func " << coef->GetName() << " is not of type RooAbsReal, fatal error" << endl ;
assert(0) ;
}
_funcList.add(*func) ;
} else {
_haveLastCoef = kTRUE ;
}
delete funcIter ;
delete coefIter ;
}
RooRealSumPdf::RooRealSumPdf(const RooRealSumPdf& other, const char* name) :
RooAbsPdf(other,name),
_normIntMgr(other._normIntMgr,this),
_haveLastCoef(other._haveLastCoef),
_funcList("!funcList",this,other._funcList),
_coefList("!coefList",this,other._coefList),
_extended(other._extended),
_doFloor(other._doFloor)
{
_funcIter = _funcList.createIterator() ;
_coefIter = _coefList.createIterator() ;
}
RooRealSumPdf::~RooRealSumPdf()
{
delete _funcIter ;
delete _coefIter ;
}
RooAbsPdf::ExtendMode RooRealSumPdf::extendMode() const
{
return (_extended && (_funcList.getSize()==_coefList.getSize())) ? CanBeExtended : CanNotBeExtended ;
}
Double_t RooRealSumPdf::evaluate() const
{
Double_t value(0) ;
RooFIter funcIter = _funcList.fwdIterator() ;
RooFIter coefIter = _coefList.fwdIterator() ;
RooAbsReal* coef ;
RooAbsReal* func ;
Double_t lastCoef(1) ;
while((coef=(RooAbsReal*)coefIter.next())) {
func = (RooAbsReal*)funcIter.next() ;
Double_t coefVal = coef->getVal() ;
if (coefVal) {
cxcoutD(Eval) << "RooRealSumPdf::eval(" << GetName() << ") coefVal = " << coefVal << " funcVal = " << func->IsA()->GetName() << "::" << func->GetName() << " = " << func->getVal() << endl ;
if (func->isSelectedComp()) {
value += func->getVal()*coefVal ;
}
lastCoef -= coef->getVal() ;
}
}
if (!_haveLastCoef) {
func = (RooAbsReal*) funcIter.next() ;
if (func->isSelectedComp()) {
value += func->getVal()*lastCoef ;
}
cxcoutD(Eval) << "RooRealSumPdf::eval(" << GetName() << ") lastCoef = " << lastCoef << " funcVal = " << func->getVal() << endl ;
if (lastCoef<0 || lastCoef>1) {
coutW(Eval) << "RooRealSumPdf::evaluate(" << GetName()
<< " WARNING: sum of FUNC coefficients not in range [0-1], value="
<< 1-lastCoef << endl ;
}
}
if (value<0 && (_doFloor || _doFloorGlobal)) {
value = 0 ;
}
return value ;
}
Bool_t RooRealSumPdf::checkObservables(const RooArgSet* nset) const
{
Bool_t ret(kFALSE) ;
_funcIter->Reset() ;
_coefIter->Reset() ;
RooAbsReal* coef ;
RooAbsReal* func ;
while((coef=(RooAbsReal*)_coefIter->Next())) {
func = (RooAbsReal*)_funcIter->Next() ;
if (func->observableOverlaps(nset,*coef)) {
coutE(InputArguments) << "RooRealSumPdf::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName()
<< " and FUNC " << func->GetName() << " have one or more observables in common" << endl ;
ret = kTRUE ;
}
if (coef->dependsOn(*nset)) {
coutE(InputArguments) << "RooRealPdf::checkObservables(" << GetName() << "): ERROR coefficient " << coef->GetName()
<< " depends on one or more of the following observables" ; nset->Print("1") ;
ret = kTRUE ;
}
}
return ret ;
}
Int_t RooRealSumPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars,
const RooArgSet* normSet2, const char* rangeName) const
{
if (allVars.getSize()==0) return 0 ;
if (_forceNumInt) return 0 ;
analVars.add(allVars) ;
RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;
Int_t sterileIdx(-1) ;
CacheElem* cache = (CacheElem*) _normIntMgr.getObj(normSet,&analVars,&sterileIdx,RooNameReg::ptr(rangeName)) ;
if (cache) {
return _normIntMgr.lastIndex()+1 ;
}
cache = new CacheElem ;
_funcIter->Reset() ;
RooAbsReal *func ;
while((func=(RooAbsReal*)_funcIter->Next())) {
RooAbsReal* funcInt = func->createIntegral(analVars,rangeName) ;
cache->_funcIntList.addOwned(*funcInt) ;
if (normSet && normSet->getSize()>0) {
RooAbsReal* funcNorm = func->createIntegral(*normSet) ;
cache->_funcNormList.addOwned(*funcNorm) ;
}
}
Int_t code = _normIntMgr.setObj(normSet,&analVars,(RooAbsCacheElement*)cache,RooNameReg::ptr(rangeName)) ;
if (normSet) {
delete normSet ;
}
return code+1 ;
}
Double_t RooRealSumPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet2, const char* rangeName) const
{
if (code==0) return getVal(normSet2) ;
CacheElem* cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ;
if (cache==0) {
std::auto_ptr<RooArgSet> vars( getParameters(RooArgSet()) );
std::auto_ptr<RooArgSet> iset( _normIntMgr.nameSet2ByIndex(code-1)->select(*vars) );
std::auto_ptr<RooArgSet> nset( _normIntMgr.nameSet1ByIndex(code-1)->select(*vars) );
RooArgSet dummy;
Int_t code2 = getAnalyticalIntegralWN(*iset,dummy,nset.get(),rangeName);
R__ASSERT(code==code2);
cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ;
R__ASSERT(cache!=0);
}
RooFIter funcIntIter = cache->_funcIntList.fwdIterator() ;
RooFIter coefIter = _coefList.fwdIterator() ;
RooFIter funcIter = _funcList.fwdIterator() ;
RooAbsReal *coef(0), *funcInt(0), *func(0) ;
Double_t value(0) ;
Double_t lastCoef(1) ;
while((coef=(RooAbsReal*)coefIter.next())) {
funcInt = (RooAbsReal*)funcIntIter.next() ;
func = (RooAbsReal*)funcIter.next() ;
Double_t coefVal = coef->getVal(normSet2) ;
if (coefVal) {
assert(func);
if (normSet2 ==0 || func->isSelectedComp()) {
assert(funcInt);
value += funcInt->getVal()*coefVal ;
}
lastCoef -= coef->getVal(normSet2) ;
}
}
if (!_haveLastCoef) {
funcInt = (RooAbsReal*) funcIntIter.next() ;
if (normSet2 ==0 || func->isSelectedComp()) {
assert(funcInt);
value += funcInt->getVal()*lastCoef ;
}
if (lastCoef<0 || lastCoef>1) {
coutW(Eval) << "RooRealSumPdf::evaluate(" << GetName()
<< " WARNING: sum of FUNC coefficients not in range [0-1], value="
<< 1-lastCoef << endl ;
}
}
Double_t normVal(1) ;
if (normSet2 && normSet2->getSize()>0) {
normVal = 0 ;
RooAbsReal* funcNorm ;
RooFIter funcNormIter = cache->_funcNormList.fwdIterator() ;
RooFIter coefIter2 = _coefList.fwdIterator() ;
while((coef=(RooAbsReal*)coefIter2.next())) {
funcNorm = (RooAbsReal*)funcNormIter.next() ;
Double_t coefVal = coef->getVal(normSet2) ;
if (coefVal) {
assert(funcNorm);
normVal += funcNorm->getVal()*coefVal ;
}
}
if (!_haveLastCoef) {
funcNorm = (RooAbsReal*) funcNormIter.next() ;
assert(funcNorm);
normVal += funcNorm->getVal()*lastCoef ;
}
}
return value / normVal;
}
Double_t RooRealSumPdf::expectedEvents(const RooArgSet* nset) const
{
Double_t n = getNorm(nset) ;
if (n<0) {
logEvalError("Expected number of events is negative") ;
}
return n ;
}
std::list<Double_t>* RooRealSumPdf::binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
list<Double_t>* sumBinB = 0 ;
Bool_t needClean(kFALSE) ;
RooFIter iter = _funcList.fwdIterator() ;
RooAbsReal* func ;
while((func=(RooAbsReal*)iter.next())) {
list<Double_t>* funcBinB = func->binBoundaries(obs,xlo,xhi) ;
if (funcBinB) {
if (!sumBinB) {
sumBinB = funcBinB ;
} else {
list<Double_t>* newSumBinB = new list<Double_t>(sumBinB->size()+funcBinB->size()) ;
merge(funcBinB->begin(),funcBinB->end(),sumBinB->begin(),sumBinB->end(),newSumBinB->begin()) ;
delete sumBinB ;
delete funcBinB ;
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 RooRealSumPdf::isBinnedDistribution(const RooArgSet& obs) const
{
RooFIter iter = _funcList.fwdIterator() ;
RooAbsReal* func ;
while((func=(RooAbsReal*)iter.next())) {
if (func->dependsOn(obs) && !func->isBinnedDistribution(obs)) {
return kFALSE ;
}
}
return kTRUE ;
}
std::list<Double_t>* RooRealSumPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
list<Double_t>* sumHint = 0 ;
Bool_t needClean(kFALSE) ;
RooFIter iter = _funcList.fwdIterator() ;
RooAbsReal* func ;
while((func=(RooAbsReal*)iter.next())) {
list<Double_t>* funcHint = func->plotSamplingHint(obs,xlo,xhi) ;
if (funcHint) {
if (!sumHint) {
sumHint = funcHint ;
} else {
list<Double_t>* newSumHint = new list<Double_t>(sumHint->size()+funcHint->size()) ;
merge(funcHint->begin(),funcHint->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 ;
}
void RooRealSumPdf::setCacheAndTrackHints(RooArgSet& trackNodes)
{
RooFIter siter = funcList().fwdIterator() ;
RooAbsArg* sarg ;
while ((sarg=siter.next())) {
if (sarg->canNodeBeCached()==Always) {
trackNodes.add(*sarg) ;
}
}
}
void RooRealSumPdf::printMetaArgs(ostream& os) const
{
_funcIter->Reset() ;
_coefIter->Reset() ;
Bool_t first(kTRUE) ;
RooAbsArg* coef, *func ;
if (_coefList.getSize()!=0) {
while((coef=(RooAbsArg*)_coefIter->Next())) {
if (!first) {
os << " + " ;
} else {
first = kFALSE ;
}
func=(RooAbsArg*)_funcIter->Next() ;
os << coef->GetName() << " * " << func->GetName() ;
}
func = (RooAbsArg*) _funcIter->Next() ;
if (func) {
os << " + [%] * " << func->GetName() ;
}
} else {
while((func=(RooAbsArg*)_funcIter->Next())) {
if (!first) {
os << " + " ;
} else {
first = kFALSE ;
}
os << func->GetName() ;
}
}
os << " " ;
}