#include "RooFit.h"
#include "TClass.h"
#include "Riostream.h"
#include "TObjString.h"
#include "TH1.h"
#include "RooRealIntegral.h"
#include "RooArgSet.h"
#include "RooAbsRealLValue.h"
#include "RooAbsCategoryLValue.h"
#include "RooRealBinding.h"
#include "RooRealAnalytic.h"
#include "RooInvTransform.h"
#include "RooSuperCategory.h"
#include "RooNumIntFactory.h"
#include "RooNumIntConfig.h"
#include "RooNameReg.h"
ClassImp(RooRealIntegral)
;
RooRealIntegral::RooRealIntegral(const char *name, const char *title,
const RooAbsReal& function, const RooArgSet& depList,
const RooArgSet* funcNormSet, const RooNumIntConfig* config,
const char* rangeName) :
RooAbsReal(name,title),
_valid(kTRUE),
_sumList("sumList","Categories to be summed numerically",this,kFALSE,kFALSE),
_intList("intList","Variables to be integrated numerically",this,kFALSE,kFALSE),
_anaList("anaList","Variables to be integrated analytically",this,kFALSE,kFALSE),
_jacList("jacList","Jacobian product term",this,kFALSE,kFALSE),
_facList("facList","Variables independent of function",this,kFALSE,kTRUE),
_facListIter(_facList.createIterator()),
_jacListIter(_jacList.createIterator()),
_function("function","Function to be integrated",this,
const_cast<RooAbsReal&>(function),kFALSE,kFALSE),
_iconfig((RooNumIntConfig*)config),
_funcACleanBranchIter(_funcACleanBranchList.createIterator()),
_mode(0),
_operMode(Hybrid),
_restartNumIntEngine(kFALSE),
_numIntEngine(0),
_numIntegrand(0),
_rangeName((TNamed*)RooNameReg::ptr(rangeName))
{
if (!_iconfig) _iconfig = (RooNumIntConfig*) function.getIntegratorConfig() ;
if (funcNormSet) {
_funcNormSet = new RooArgSet ;
TIterator* iter = funcNormSet->createIterator() ;
RooAbsArg* nArg ;
while ((nArg=(RooAbsArg*)iter->Next())) {
if (function.dependsOn(*nArg)) {
_funcNormSet->addClone(*nArg) ;
}
}
delete iter ;
} else {
_funcNormSet = 0 ;
}
RooArgSet intDepList(depList) ;
RooAbsArg *arg ;
TIterator* depIter = intDepList.createIterator() ;
while((arg=(RooAbsArg*)depIter->Next())) {
if(!arg->isLValue()) {
cout << ClassName() << "::" << GetName() << ": cannot integrate non-lvalue ";
arg->Print("1");
_valid= kFALSE;
}
if (!function.dependsOn(*arg)) {
RooAbsArg* argClone = (RooAbsArg*) arg->Clone() ;
_facListOwned.addOwned(*argClone) ;
_facList.add(*argClone) ;
addServer(*argClone,kFALSE,kTRUE) ;
}
}
RooArgSet exclLVBranches("exclLVBranches") ;
RooArgSet branchList ;
function.branchNodeServerList(&branchList) ;
TIterator* bIter = branchList.createIterator() ;
RooAbsArg* branch ;
while((branch=(RooAbsArg*)bIter->Next())) {
RooAbsRealLValue *realArgLV = dynamic_cast<RooAbsRealLValue*>(branch) ;
RooAbsCategoryLValue *catArgLV = dynamic_cast<RooAbsCategoryLValue*>(branch) ;
if ((realArgLV && (realArgLV->isJacobianOK(intDepList)!=0)) || catArgLV) {
exclLVBranches.add(*branch) ;
}
}
delete bIter ;
exclLVBranches.remove(depList,kTRUE,kTRUE) ;
RooArgSet exclLVServers("exclLVServers") ;
exclLVServers.add(intDepList) ;
TIterator *sIter = exclLVServers.createIterator() ;
bIter = exclLVBranches.createIterator() ;
RooAbsArg *server ;
Bool_t converged(kFALSE) ;
while(!converged) {
converged=kTRUE ;
sIter->Reset() ;
while ((server=(RooAbsArg*)sIter->Next())) {
if (!servesExclusively(server,exclLVBranches)) {
exclLVServers.remove(*server) ;
converged=kFALSE ;
}
}
bIter->Reset() ;
while((branch=(RooAbsArg*)bIter->Next())) {
RooArgSet* brDepList = branch->getObservables(&intDepList) ;
RooArgSet bsList(*brDepList,"bsList") ;
delete brDepList ;
bsList.remove(exclLVServers,kTRUE,kTRUE) ;
if (bsList.getSize()>0) {
exclLVBranches.remove(*branch,kTRUE,kTRUE) ;
converged=kFALSE ;
}
}
}
delete sIter ;
delete bIter ;
if (exclLVServers.getSize()>0) {
intDepList.remove(exclLVServers) ;
intDepList.add(exclLVBranches) ;
}
RooArgSet anIntOKDepList ;
depIter->Reset() ;
while((arg=(RooAbsArg*)depIter->Next())) {
if (function.forceAnalyticalInt(*arg)) {
anIntOKDepList.add(*arg) ;
}
}
delete depIter ;
sIter = function.serverIterator() ;
while((arg=(RooAbsArg*)sIter->Next())) {
if (!arg->dependsOn(intDepList)) {
addServer(*arg,kTRUE,kFALSE) ;
continue ;
} else {
RooArgSet argLeafServers ;
arg->leafNodeServerList(&argLeafServers) ;
TIterator* lIter = argLeafServers.createIterator() ;
RooAbsArg* leaf ;
while((leaf=(RooAbsArg*)lIter->Next())) {
if (depList.find(leaf->GetName())) {
addServer(*leaf,kFALSE,kTRUE) ;
} else {
addServer(*leaf,kTRUE,kFALSE) ;
}
}
delete lIter ;
}
Bool_t depOK(kFALSE) ;
if (arg->isDerived()) {
RooAbsRealLValue *realArgLV = dynamic_cast<RooAbsRealLValue*>(arg) ;
RooAbsCategoryLValue *catArgLV = dynamic_cast<RooAbsCategoryLValue*>(arg) ;
if ((realArgLV && intDepList.find(realArgLV->GetName()) && (realArgLV->isJacobianOK(intDepList)!=0)) || catArgLV) {
depOK = kTRUE ;
Bool_t overlapOK = kTRUE ;
RooAbsArg *otherArg ;
TIterator* sIter2 = function.serverIterator() ;
while((otherArg=(RooAbsArg*)sIter2->Next())) {
if (arg==otherArg) continue ;
if (arg->overlaps(*otherArg)) {
overlapOK=kFALSE ;
}
}
if (!overlapOK) depOK=kFALSE ;
delete sIter2 ;
}
} else {
depOK = kTRUE ;
}
if (depOK) {
anIntOKDepList.add(*arg,kTRUE) ;
}
}
RooArgSet anIntDepList ;
_mode = ((RooAbsReal&)_function.arg()).getAnalyticalIntegralWN(anIntOKDepList,_anaList,_funcNormSet,RooNameReg::str(_rangeName)) ;
if (_mode==0) {
_anaList.removeAll() ;
}
function.getVal(_funcNormSet) ;
RooArgSet numIntDepList ;
TIterator* aiIter = _anaList.createIterator() ;
while ((arg=(RooAbsArg*)aiIter->Next())) {
if (arg->IsA()->InheritsFrom(RooAbsRealLValue::Class()) && arg->isDerived()) {
_jacList.add(*arg) ;
RooAbsArg *argDep ;
RooArgSet *argDepList = arg->getObservables(&intDepList) ;
TIterator *adIter = argDepList->createIterator() ;
while ((argDep=(RooAbsArg*)adIter->Next())) {
if (argDep->IsA()->InheritsFrom(RooAbsCategoryLValue::Class()) && intDepList.contains(*argDep)) {
numIntDepList.add(*argDep,kTRUE) ;
}
}
delete adIter ;
delete argDepList ;
}
}
delete aiIter ;
sIter->Reset() ;
while((arg=(RooAbsArg*)sIter->Next())) {
if (!_anaList.find(arg->GetName()) && arg->dependsOn(intDepList)) {
if (dynamic_cast<RooAbsLValue*>(arg) && arg->isDerived() && intDepList.contains(*arg)) {
numIntDepList.add(*arg,kTRUE) ;
} else {
RooArgSet *argDeps = arg->getObservables(&intDepList) ;
TIterator* iter = argDeps->createIterator() ;
RooAbsArg* dep ;
while((dep=(RooAbsArg*)iter->Next())) {
if (!_anaList.find(dep->GetName())) {
numIntDepList.add(*dep,kTRUE) ;
}
}
delete iter ;
delete argDeps ;
}
}
}
delete sIter ;
TIterator* numIter=numIntDepList.createIterator() ;
while ((arg=(RooAbsArg*)numIter->Next())) {
if (arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
_intList.add(*arg) ;
} else if (arg->IsA()->InheritsFrom(RooAbsCategoryLValue::Class())) {
_sumList.add(*arg) ;
}
}
delete numIter ;
if (numIntDepList.getSize()>0) {
_operMode = Hybrid ;
} else if (_anaList.getSize()>0) {
_operMode = Analytic ;
} else {
_operMode = PassThrough ;
}
autoSelectDirtyMode() ;
}
void RooRealIntegral::autoSelectDirtyMode()
{
TIterator* siter = serverIterator() ;
RooAbsArg* server ;
while((server=(RooAbsArg*)siter->Next())){
RooArgSet leafSet ;
server->leafNodeServerList(&leafSet) ;
TIterator* liter = leafSet.createIterator() ;
RooAbsArg* leaf ;
while((leaf=(RooAbsArg*)liter->Next())) {
if (leaf->operMode()==ADirty && leaf->isValueServer(*this)) {
setOperMode(ADirty) ;
break ;
}
if (leaf->getAttribute("projectedDependent")) {
setOperMode(ADirty) ;
break ;
}
}
delete liter ;
}
delete siter ;
}
Bool_t RooRealIntegral::servesExclusively(const RooAbsArg* server,const RooArgSet& exclLVBranches) const
{
if (exclLVBranches.getSize()==0) return kFALSE ;
if (server->_clientList.GetSize()==0 && exclLVBranches.find(server->GetName())) {
return kFALSE ;
}
Int_t numLVServ(0) ;
RooAbsArg* client ;
TIterator* cIter = server->clientIterator() ;
while((client=(RooAbsArg*)cIter->Next())) {
if (!exclLVBranches.find(client->GetName())) {
if (!servesExclusively(client,exclLVBranches)) {
return kFALSE ;
}
} else {
numLVServ++ ;
}
}
delete cIter ;
return (numLVServ==1) ;
}
Bool_t RooRealIntegral::initNumIntegrator() const
{
if(0 != _numIntEngine) {
if(_numIntEngine->isValid() && _numIntEngine->checkLimits() && !_restartNumIntEngine ) return kTRUE;
delete _numIntEngine ;
_numIntEngine= 0;
if(0 != _numIntegrand) {
delete _numIntegrand;
_numIntegrand= 0;
}
}
if(0 == _intList.getSize()) return kTRUE;
if(_mode != 0) {
_numIntegrand= new RooRealAnalytic(_function.arg(),_intList,_mode,_funcNormSet,_rangeName);
}
else {
_numIntegrand= new RooRealBinding(_function.arg(),_intList,_funcNormSet,kFALSE,_rangeName);
}
if(0 == _numIntegrand || !_numIntegrand->isValid()) {
cout << ClassName() << "::" << GetName() << ": failed to create valid integrand." << endl;
return kFALSE;
}
_numIntEngine = RooNumIntFactory::instance().createIntegrator(*_numIntegrand,*_iconfig) ;
if(0 == _numIntEngine || !_numIntEngine->isValid()) {
cout << ClassName() << "::" << GetName() << ": failed to create valid integrator." << endl;
return kFALSE;
}
_restartNumIntEngine = kFALSE ;
return kTRUE;
}
RooRealIntegral::RooRealIntegral(const RooRealIntegral& other, const char* name) :
RooAbsReal(other,name),
_valid(other._valid),
_sumList("sumList",this,other._sumList),
_intList("intList",this,other._intList),
_anaList("anaList",this,other._anaList),
_jacList("jacList",this,other._jacList),
_facList("facList","Variables independent of function",this,kFALSE,kTRUE),
_facListIter(_facList.createIterator()),
_jacListIter(_jacList.createIterator()),
_function("function",this,other._function),
_iconfig(other._iconfig),
_funcACleanBranchIter(_funcACleanBranchList.createIterator()),
_mode(other._mode),
_operMode(other._operMode),
_restartNumIntEngine(kFALSE),
_numIntEngine(0),
_numIntegrand(0),
_rangeName(other._rangeName)
{
_funcNormSet = other._funcNormSet ? (RooArgSet*)other._funcNormSet->snapshot(kFALSE) : 0 ;
other._facListIter->Reset() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)other._facListIter->Next())) {
RooAbsArg* argClone = (RooAbsArg*) arg->Clone() ;
_facListOwned.addOwned(*argClone) ;
_facList.add(*argClone) ;
addServer(*argClone,kFALSE,kTRUE) ;
}
}
RooRealIntegral::~RooRealIntegral()
{
if (_numIntEngine) delete _numIntEngine ;
if (_numIntegrand) delete _numIntegrand ;
if (_funcNormSet) delete _funcNormSet ;
delete _funcACleanBranchIter ;
delete _facListIter ;
delete _jacListIter ;
}
Double_t RooRealIntegral::evaluate() const
{
Double_t retVal(0) ;
switch (_operMode) {
case Hybrid:
{
prepareACleanFunc() ;
if(!(_valid= initNumIntegrator())) {
cout << ClassName() << "::" << GetName()
<< ":evaluate: cannot initialize numerical integrator" << endl;
return 0;
}
RooArgSet *saveInt = (RooArgSet*) _intList.snapshot() ;
RooArgSet *saveSum = (RooArgSet*) _sumList.snapshot() ;
retVal = sum() / jacobianProduct() ;
_intList=*saveInt ;
_sumList=*saveSum ;
delete saveInt ;
delete saveSum ;
restoreACleanFunc() ;
break ;
}
case Analytic:
{
retVal = ((RooAbsReal&)_function.arg()).analyticalIntegralWN(_mode,_funcNormSet,RooNameReg::str(_rangeName)) / jacobianProduct() ;
if (RooAbsPdf::_verboseEval>0)
cout << "RooRealIntegral::evaluate_analytic(" << GetName()
<< ")func = " << _function.arg().IsA()->GetName() << "::" << _function.arg().GetName()
<< " raw = " << retVal << endl ;
break ;
}
case PassThrough:
{
retVal= _function.arg().getVal(_funcNormSet) ;
break ;
}
}
RooAbsArg *arg ;
_facListIter->Reset() ;
while((arg=(RooAbsArg*)_facListIter->Next())) {
if (arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
RooAbsRealLValue* argLV = (RooAbsRealLValue*)arg ;
retVal *= (argLV->getMax() - argLV->getMin()) ;
}
if (arg->IsA()->InheritsFrom(RooAbsCategoryLValue::Class())) {
RooAbsCategoryLValue* argLV = (RooAbsCategoryLValue*)arg ;
retVal *= argLV->numTypes() ;
}
}
if (RooAbsPdf::_verboseEval>0) {
cout << "RooRealIntegral::evaluate(" << GetName() << ") raw*fact = " << retVal << endl ;
}
return retVal ;
}
void RooRealIntegral::prepareACleanFunc() const
{
if (_funcBranchList.getSize()==0) {
_function.arg().branchNodeServerList(&_funcBranchList) ;
}
_funcACleanBranchList.removeAll() ;
_funcACleanBranchList.add(_funcBranchList) ;
_funcACleanBranchIter->Reset() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)_funcACleanBranchIter->Next())) {
if (arg->operMode()!=RooAbsArg::AClean) {
_funcACleanBranchList.remove(*arg) ;
} else {
arg->setOperMode(RooAbsArg::ADirty, kFALSE) ;
}
}
}
void RooRealIntegral::restoreACleanFunc() const
{
_funcACleanBranchIter->Reset() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)_funcACleanBranchIter->Next())) {
arg->setOperMode(RooAbsArg::AClean) ;
}
}
Double_t RooRealIntegral::jacobianProduct() const
{
Double_t jacProd(1) ;
_jacListIter->Reset() ;
RooAbsRealLValue* arg ;
while ((arg=(RooAbsRealLValue*)_jacListIter->Next())) {
jacProd *= arg->jacobian() ;
}
return jacProd ;
}
Double_t RooRealIntegral::sum() const
{
if (_sumList.getSize()!=0) {
Double_t total(0) ;
RooSuperCategory sumCat("sumCat","sumCat",_sumList) ;
TIterator* sumIter = sumCat.typeIterator() ;
RooCatType* type ;
while((type=(RooCatType*)sumIter->Next())) {
sumCat.setIndex(type->getVal()) ;
if (!_rangeName || sumCat.inRange(RooNameReg::str(_rangeName))) {
total += integrate() / jacobianProduct() ;
}
}
delete sumIter ;
return total ;
} else {
Double_t ret = integrate() / jacobianProduct() ;
return ret ;
}
}
Double_t RooRealIntegral::integrate() const
{
if (!_numIntEngine) {
return ((RooAbsReal&)_function.arg()).analyticalIntegralWN(_mode,_funcNormSet,RooNameReg::str(_rangeName)) ;
}
else {
return _numIntEngine->calculate() ;
}
}
Bool_t RooRealIntegral::redirectServersHook(const RooAbsCollection& ,
Bool_t , Bool_t , Bool_t )
{
_funcBranchList.removeAll() ;
_restartNumIntEngine = kTRUE ;
autoSelectDirtyMode() ;
return kFALSE ;
}
Bool_t RooRealIntegral::isValidReal(Double_t , Bool_t ) const
{
return kTRUE ;
}
void RooRealIntegral::printToStream(ostream& os, PrintOption opt, TString indent) const
{
RooAbsReal::printToStream(os,opt,indent) ;
if (opt==Verbose) {
os << indent << "--- RooRealIntegral ---" << endl;
os << indent << " Integrates ";
_function.arg().printToStream(os,Standard,indent);
TString deeper(indent);
deeper.Append(" ");
os << indent << " operating mode is "
<< (_operMode==Hybrid?"Hybrid":(_operMode==Analytic?"Analytic":"PassThrough")) << endl ;
os << indent << " Summed discrete args are ";
_sumList.printToStream(os,Standard,deeper);
os << indent << " Numerically integrated args are ";
_intList.printToStream(os,Standard,deeper);
os << indent << " Analytically integrated args using mode " << _mode << " are ";
_anaList.printToStream(os,Standard,deeper);
os << indent << " Arguments included in Jacobian are ";
_jacList.printToStream(os,Standard,deeper);
os << indent << " Factorized arguments are ";
_facList.printToStream(os,Standard,deeper);
os << indent << " Function normalization set " ;
if (_funcNormSet)
_funcNormSet->Print("1") ;
else
os << "<none>" ;
os << endl ;
return ;
}
}
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.