RooRealIntegral.cxx

Go to the documentation of this file.

/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * 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)             *
 *****************************************************************************/

/**
\file RooRealIntegral.cxx
\class RooRealIntegral
\ingroup Roofitcore

RooRealIntegral performs hybrid numerical/analytical integrals of RooAbsReal objects
The class performs none of the actual integration, but only manages the logic
of what variables can be integrated analytically, accounts for eventual jacobian
terms and defines what numerical integrations needs to be done to complement the
analytical integral.
The actual analytical integrations (if any) are done in the PDF themselves, the numerical
integration is performed in the various implemenations of the RooAbsIntegrator base class.
**/

#include "RooFit.h"

#include "TClass.h"
#include "RooMsgService.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"
#include "RooExpensiveObjectCache.h"
#include "RooConstVar.h"
#include "RooDouble.h"
#include "RooTrace.h"

using namespace std;

ClassImp(RooRealIntegral); 
;


Int_t RooRealIntegral::_cacheAllNDim(2) ;


////////////////////////////////////////////////////////////////////////////////

RooRealIntegral::RooRealIntegral() : 
  _valid(kFALSE),
  _respectCompSelect(kFALSE),
  _funcNormSet(0),
  _iconfig(0),
  _sumCatIter(0),
  _mode(0),
  _intOperMode(Hybrid),
  _restartNumIntEngine(kFALSE),
  _numIntEngine(0),
  _numIntegrand(0),
  _rangeName(0),
  _params(0),
  _cacheNum(kFALSE)
{
  _facListIter = _facList.createIterator() ;
  _jacListIter = _jacList.createIterator() ;
  TRACE_CREATE
}



////////////////////////////////////////////////////////////////////////////////
/// Construct integral of 'function' over observables in 'depList'
/// in range 'rangeName'  with normalization observables 'funcNormSet' 
/// (for p.d.f.s). In the integral is performed to the maximum extent
/// possible the internal (analytical) integrals advertised by function.
/// The other integrations are performed numerically. The optional
/// config object prescribes how these numeric integrations are configured.
///

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), 
  _respectCompSelect(kFALSE),
  _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("!func","Function to be integrated",this,
       const_cast<RooAbsReal&>(function),kFALSE,kFALSE), 
  _iconfig((RooNumIntConfig*)config),
  _sumCat("!sumCat","SuperCategory for summation",this,kFALSE,kFALSE),
  _sumCatIter(0),
  _mode(0),
  _intOperMode(Hybrid), 
  _restartNumIntEngine(kFALSE),
  _numIntEngine(0), 
  _numIntegrand(0),
  _rangeName((TNamed*)RooNameReg::ptr(rangeName)),
  _params(0),
  _cacheNum(kFALSE)
{
  //   A) Check that all dependents are lvalues 
  //
  //   B) Check if list of dependents can be re-expressed in        
  //      lvalues that are higher in the expression tree            
  //
  //   C) Check for dependents that the PDF insists on integrating  
  //      analytically iself                                        
  //
  //   D) Make list of servers that can be integrated analytically  
  //      Add all parameters/dependents as value/shape servers      
  //
  //   E) Interact with function to make list of objects actually integrated analytically   
  //
  //   F) Make list of numerical integration variables consisting of:               
  //     - Category dependents of RealLValues in analytical integration             
  //     - Leaf nodes server lists of function server that are not analytically integrated   
  //     - Make Jacobian list for analytically integrated RealLValues            
  //
  //   G) Split numeric list in integration list and summation list   
  //

  oocxcoutI(&function,Integration) << "RooRealIntegral::ctor(" << GetName() << ") Constructing integral of function " 
                 << function.GetName() << " over observables" << depList << " with normalization " 
                 << (funcNormSet?*funcNormSet:RooArgSet()) << " with range identifier " 
                 << (rangeName?rangeName:"<none>") << endl ;

  
  // Choose same expensive object cache as integrand
  setExpensiveObjectCache(function.expensiveObjectCache()) ;
//   cout << "RRI::ctor(" << GetName() << ") setting expensive object cache to " << &expensiveObjectCache() << " as taken from " << function.GetName() << endl ;

  // Use objects integrator configuration if none is specified
  if (!_iconfig) _iconfig = (RooNumIntConfig*) function.getIntegratorConfig() ;

  // Save private copy of funcNormSet, if supplied, excluding factorizing terms
  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 ;
  }

  //_funcNormSet = funcNormSet ? (RooArgSet*)funcNormSet->snapshot(kFALSE) : 0 ;
  
  // Make internal copy of dependent list
  RooArgSet intDepList(depList) ;

  RooAbsArg *arg ;

  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * A) Check that all dependents are lvalues and filter out any
  //      dependents that the PDF doesn't explicitly depend on
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  
  TIterator* depIter = intDepList.createIterator() ;
  while((arg=(RooAbsArg*)depIter->Next())) {
    if(!arg->isLValue()) {
      coutE(InputArguments) << 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) ;
    }
  }

  if (_facList.getSize()>0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Factorizing obserables are " << _facList << endl ;
  }


  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * B) Check if list of dependents can be re-expressed in       *
  // *    lvalues that are higher in the expression tree           *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 


  // Initial fill of list of LValue branches
  RooArgSet exclLVBranches("exclLVBranches") ;
  RooArgSet branchList,branchListVD ;
  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) ;
//       cout << "exclv branch = " << endl ;
//       branch->printCompactTree() ;
    }
    if (dependsOnValue(*branch)) {
      branchListVD.add(*branch) ;
    } else {
//       cout << "value of self does not depend on branch " << branch->GetName() << endl ;
    }
  }
  delete bIter ;
  exclLVBranches.remove(depList,kTRUE,kTRUE) ;
//    cout << "exclLVBranches = " << exclLVBranches << endl ;

  // Initial fill of list of LValue leaf servers (put in intDepList)
  RooArgSet exclLVServers("exclLVServers") ;
  exclLVServers.add(intDepList) ;

//    cout << "begin exclLVServers = " << exclLVServers << endl ;
  
  // Obtain mutual exclusive dependence by iterative reduction
  TIterator *sIter = exclLVServers.createIterator() ;
  bIter = exclLVBranches.createIterator() ;
  RooAbsArg *server ;
  Bool_t converged(kFALSE) ;
  while(!converged) {
    converged=kTRUE ;

    // Reduce exclLVServers to only those serving exclusively exclLVBranches
    sIter->Reset() ;
    while ((server=(RooAbsArg*)sIter->Next())) {
      if (!servesExclusively(server,exclLVBranches,branchListVD)) {
   exclLVServers.remove(*server) ;
//    cout << "removing " << server->GetName() << " from exclLVServers because servesExclusively(" << server->GetName() << "," << exclLVBranches << ") faile" << endl ;
   converged=kFALSE ;
      }
    }
    
    // Reduce exclLVBranches to only those depending exclusisvely on exclLVservers
    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) ;
//    cout << "removing " << branch->GetName() << " from exclLVBranches" << endl ;
   converged=kFALSE ;
      }
    }
  }

  // Eliminate exclLVBranches that do not depend on any LVServer
  bIter->Reset() ;
  while((branch=(RooAbsArg*)bIter->Next())) {    
    if (!branch->dependsOnValue(exclLVServers)) {
      //cout << "LV branch " << branch->GetName() << " does not depend on any LVServer (" << exclLVServers << ") and will be removed" << endl ; 
      exclLVBranches.remove(*branch,kTRUE,kTRUE) ;
    }
  } 

  delete sIter ;
  delete bIter ;

//   cout << "end exclLVServers = " << exclLVServers << endl ;
     
  // Replace exclusive lvalue branch servers with lvalue branches
  // WVE Don't do this for binned distributions - deal with this using numeric integration with transformed bin boundaroes
  if (exclLVServers.getSize()>0 && !function.isBinnedDistribution(exclLVBranches)) {
//     cout << "activating LVservers " << exclLVServers << " for use in integration " << endl ;
    intDepList.remove(exclLVServers) ;
    intDepList.add(exclLVBranches) ;

    //cout << "intDepList removing exclLVServers " << exclLVServers << endl ;
    //cout << "intDepList adding exclLVBranches " << exclLVBranches << endl ;

  }

     
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * C) Check for dependents that the PDF insists on integrating *
  //      analytically iself                                       *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  RooArgSet anIntOKDepList ;
  depIter->Reset() ;
  while((arg=(RooAbsArg*)depIter->Next())) {
    if (function.forceAnalyticalInt(*arg)) {
      anIntOKDepList.add(*arg) ;
    }
  }
  delete depIter ;
  
  if (anIntOKDepList.getSize()>0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Observables that function forcibly requires to be integrated internally " << anIntOKDepList << endl ;
  }


  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * D) Make list of servers that can be integrated analytically *
  //      Add all parameters/dependents as value/shape servers     *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  sIter = function.serverIterator() ;
  while((arg=(RooAbsArg*)sIter->Next())) {

    //cout << "considering server" << arg->GetName() << endl ;
    
    // Dependent or parameter?
    if (!arg->dependsOnValue(intDepList)) {

      //cout << " server does not depend on observables, adding server as value server to integral" << endl ;

      if (function.dependsOnValue(*arg)) {
   addServer(*arg,kTRUE,kFALSE) ;
      }

      continue ;

    } else {

      // Add final dependents of arg as shape servers
      RooArgSet argLeafServers ;
      arg->leafNodeServerList(&argLeafServers,0,kFALSE) ;

      //arg->printCompactTree() ;
      //cout << "leaf nodes of server are " << argLeafServers << " depList = " << depList << endl ;

      // Skip arg if it is neither value or shape server
      if (!arg->isValueServer(function) && !arg->isShapeServer(function)) {
   //cout << " server is neither value not shape server of function, ignoring" << endl ;
   continue ;
      }
      
      TIterator* lIter = argLeafServers.createIterator() ;
      RooAbsArg* leaf ;
      while((leaf=(RooAbsArg*)lIter->Next())) {

   //cout << " considering leafnode " << leaf->GetName() << " of server " << arg->GetName() << endl ;

   if (depList.find(leaf->GetName()) && function.dependsOnValue(*leaf)) {

     RooAbsRealLValue* leaflv = dynamic_cast<RooAbsRealLValue*>(leaf) ;
     if (leaflv && leaflv->getBinning(rangeName).isParameterized()) {
       oocxcoutD(&function,Integration) << function.GetName() << " : Observable " << leaf->GetName() << " has parameterized binning, add value dependence of boundary objects rather than shape of leaf" << endl ;
       if (leaflv->getBinning(rangeName).lowBoundFunc()) {
         addServer(*leaflv->getBinning(rangeName).lowBoundFunc(),kTRUE,kFALSE) ;
       }
       if(leaflv->getBinning(rangeName).highBoundFunc()) {
         addServer(*leaflv->getBinning(rangeName).highBoundFunc(),kTRUE,kFALSE) ;
       }
     } else {
       oocxcoutD(&function,Integration) << function.GetName() << ": Adding observable " << leaf->GetName() << " of server " 
                    << arg->GetName() << " as shape dependent" << endl ;
       addServer(*leaf,kFALSE,kTRUE) ;
     }
   } else if (!depList.find(leaf->GetName())) {

     if (function.dependsOnValue(*leaf)) {
       oocxcoutD(&function,Integration) << function.GetName() << ": Adding parameter " << leaf->GetName() << " of server " << arg->GetName() << " as value dependent" << endl ;
       addServer(*leaf,kTRUE,kFALSE) ;
     } else {
       oocxcoutD(&function,Integration) << function.GetName() << ": Adding parameter " << leaf->GetName() << " of server " << arg->GetName() << " as shape dependent" << endl ;
       addServer(*leaf,kFALSE,kTRUE) ;
     }
   }  
      }
      delete lIter ;
    }

    // If this dependent arg is self-normalized, stop here
    //if (function.selfNormalized()) continue ;

    Bool_t depOK(kFALSE) ;
    // Check for integratable AbsRealLValue

    //cout << "checking server " << arg->IsA()->GetName() << "::" << arg->GetName() << endl ;

    if (arg->isDerived()) {
      RooAbsRealLValue    *realArgLV = dynamic_cast<RooAbsRealLValue*>(arg) ;
      RooAbsCategoryLValue *catArgLV = dynamic_cast<RooAbsCategoryLValue*>(arg) ;
      //cout << "realArgLV = " << realArgLV << " intDepList = " << intDepList << endl ;
      if ((realArgLV && intDepList.find(realArgLV->GetName()) && (realArgLV->isJacobianOK(intDepList)!=0)) || catArgLV) {  

   //cout  << " arg " << arg->GetName() << " is derived LValue with valid jacobian" << endl ;

   // Derived LValue with valid jacobian
   depOK = kTRUE ;
   
   // Now, check for overlaps
   Bool_t overlapOK = kTRUE ;
   RooAbsArg *otherArg ;
   TIterator* sIter2 = function.serverIterator() ; 
   while((otherArg=(RooAbsArg*)sIter2->Next())) {
     // skip comparison with self
     if (arg==otherArg) continue ;
     if (otherArg->IsA()==RooConstVar::Class()) continue ;
     if (arg->overlaps(*otherArg,kTRUE)) {
       //cout << "arg " << arg->GetName() << " overlaps with " << otherArg->GetName() << endl ;
       //overlapOK=kFALSE ;
     }
   }        
   // coverity[DEADCODE]
   if (!overlapOK) depOK=kFALSE ;      
   
   //cout << "overlap check returns OK=" << (depOK?"T":"F") << endl ;

   delete sIter2 ;
      }
    } else {
      // Fundamental types are always OK
      depOK = kTRUE ;
    }
    
    // Add server to list of dependents that are OK for analytical integration
    if (depOK) {
      anIntOKDepList.add(*arg,kTRUE) ;      
      oocxcoutI(&function,Integration) << function.GetName() << ": Observable " << arg->GetName() << " is suitable for analytical integration (if supported by p.d.f)" << endl ;
    }
  }
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * E) interact with function to make list of objects actually integrated analytically  *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  RooArgSet anIntDepList ;

  RooArgSet *anaSet = new RooArgSet( _anaList, Form("UniqueCloneOf_%s",_anaList.GetName()));
  _mode = ((RooAbsReal&)_function.arg()).getAnalyticalIntegralWN(anIntOKDepList,*anaSet,_funcNormSet,RooNameReg::str(_rangeName)) ;    
  _anaList.removeAll() ;
  _anaList.add(*anaSet);    
  delete anaSet;

  // Avoid confusion -- if mode is zero no analytical integral is defined regardless of contents of _anaListx
  if (_mode==0) {
    _anaList.removeAll() ;
  }

  if (_mode!=0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Function integrated observables " << _anaList << " internally with code " << _mode << endl ;
  }


  // WVE kludge: synchronize dset for use in analyticalIntegral
  function.getVal(_funcNormSet) ;

  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * F) Make list of numerical integration variables consisting of:            *  
  // *   - Category dependents of RealLValues in analytical integration          *  
  // *   - Expanded server lists of server that are not analytically integrated  *
  // *    Make Jacobian list with analytically integrated RealLValues            *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  RooArgSet numIntDepList ;


  // Loop over actually analytically integrated dependents
  TIterator* aiIter = _anaList.createIterator() ;
  while ((arg=(RooAbsArg*)aiIter->Next())) {    

    // Process only derived RealLValues
    if (arg->IsA()->InheritsFrom(RooAbsRealLValue::Class()) && arg->isDerived() && !arg->isFundamental()) {

      // Add to list of Jacobians to calculate
      _jacList.add(*arg) ;

      // Add category dependent of LValueReal used in integration
      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 ;

  
  // If nothing was integrated analytically, swap back LVbranches for LVservers for subsequent numeric integration
  if (_anaList.getSize()==0) {    
    if (exclLVServers.getSize()>0) {
      //cout << "NUMINT phase analList is empty. exclLVServers = " << exclLVServers << endl ;
      intDepList.remove(exclLVBranches) ;      
      intDepList.add(exclLVServers) ;
     }            
  }
  //cout << "NUMINT intDepList = " << intDepList << endl ;

  // Loop again over function servers to add remaining numeric integrations
  sIter->Reset() ;
  while((arg=(RooAbsArg*)sIter->Next())) {

    //cout << "processing server for numeric integration " << arg->IsA()->GetName() << "::" << arg->GetName() << endl ;

    // Process only servers that are not treated analytically
    if (!_anaList.find(arg->GetName()) && arg->dependsOn(intDepList)) {

      // Process only derived RealLValues
      if (dynamic_cast<RooAbsLValue*>(arg) && arg->isDerived() && intDepList.contains(*arg)) {
   numIntDepList.add(*arg,kTRUE) ;  
      } else {

   // WVE this will only get the observables, but not l-value transformations
   // Expand server in final dependents 
   RooArgSet *argDeps = arg->getObservables(&intDepList) ;

   if (argDeps->getSize()>0) {

     // Add final dependents, that are not forcibly integrated analytically, 
     // to numerical integration list      
     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 ;

  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  // * G) Split numeric list in integration list and summation list  *
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

  // Split numeric integration list in summation and integration lists
  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 (_anaList.getSize()>0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Observables " << _anaList << " are analytically integrated with code " << _mode << endl ;
  }
  if (_intList.getSize()>0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Observables " << _intList << " are numerically integrated" << endl ;
  }
  if (_sumList.getSize()>0) {
    oocxcoutI(&function,Integration) << function.GetName() << ": Observables " << _sumList << " are numerically summed" << endl ;
  }
  

  // Determine operating mode
  if (numIntDepList.getSize()>0) {
    // Numerical and optional Analytical integration
    _intOperMode = Hybrid ;
  } else if (_anaList.getSize()>0) {
    // Purely analytical integration
    _intOperMode = Analytic ;    
  } else {
    // No integration performed
    _intOperMode = PassThrough ;
  }

  // Determine auto-dirty status  
  autoSelectDirtyMode() ;

  // Create value caches for _intList and _sumList
  _intList.snapshot(_saveInt) ;
  _sumList.snapshot(_saveSum) ;

  
  if (_sumList.getSize()>0) {
    RooSuperCategory *sumCat = new RooSuperCategory(Form("%s_sumCat",GetName()),"sumCat",_sumList) ;
    _sumCatIter = sumCat->typeIterator() ;    
    _sumCat.addOwned(*sumCat) ;
  }

  TRACE_CREATE
}



////////////////////////////////////////////////////////////////////////////////
/// Set appropriate cache operation mode for integral depending on cache operation
/// mode of server objects

void RooRealIntegral::autoSelectDirtyMode() 
{
  // If any of our servers are is forcedDirty or a projectedDependent, then we need to be ADirty
  TIterator* siter = serverIterator() ;  
  RooAbsArg* server ;
  while((server=(RooAbsArg*)siter->Next())){
    if (server->isValueServer(*this)) {
      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 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function that returns true if 'object server' is a server
/// to exactly one of the RooAbsArgs in 'exclLVBranches'

Bool_t RooRealIntegral::servesExclusively(const RooAbsArg* server,const RooArgSet& exclLVBranches, const RooArgSet& allBranches) const
{  
  // Determine if given server serves exclusively exactly one of the given nodes in exclLVBranches

  // Special case, no LV servers available
  if (exclLVBranches.getSize()==0) return kFALSE ;

  // If server has no clients and is not an LValue itself, return false
   if (server->_clientList.GetSize()==0 && exclLVBranches.find(server->GetName())) {
     return kFALSE ;
   }

   // WVE must check for value relations only here!!!!


//    cout << "servesExclusively: does " << server->GetName() << " serve only one of " << exclLVBranches << endl ;

   // Loop over all clients
   Int_t numLVServ(0) ;
   RooAbsArg* client ;
   TIterator* cIter = server->valueClientIterator() ;
   while((client=(RooAbsArg*)cIter->Next())) {
//      cout << "now checking value client " << client->GetName() << " of server " << server->GetName() << endl ;
     // If client is not an LValue, recurse
     if (!(exclLVBranches.find(client->GetName())==client)) {
//        cout << " client " << client->GetName() << "is not an lvalue" << endl ;
       if (allBranches.find(client->GetName())==client) {
//     cout << " ... recursing call" << endl ;
    if (!servesExclusively(client,exclLVBranches,allBranches)) {
    // Client is a non-LValue that doesn't have an exclusive LValue server
    delete cIter ;
//     cout << "client " << client->GetName() << " is a non-lvalue that doesn't have an exclusive lvalue server" << endl ;
    return kFALSE ;   
    }
       }
     } else {
       // Client is an LValue       
//        cout << "client " << client->GetName() << " of server " << server->GetName() << " is an LValue " << endl ;
       numLVServ++ ;
     }
   }

   delete cIter ;
//    cout << "numLVserv = " << numLVServ << endl ;
   return (numLVServ==1) ;
}




////////////////////////////////////////////////////////////////////////////////
/// (Re)Initialize numerical integration engine if necessary. Return kTRUE if
/// successful, or otherwise kFALSE.

Bool_t RooRealIntegral::initNumIntegrator() const
{
  // if we already have an engine, check if it still works for the present limits.
  if(0 != _numIntEngine) {
    if(_numIntEngine->isValid() && _numIntEngine->checkLimits() && !_restartNumIntEngine ) return kTRUE;
    // otherwise, cleanup the old engine
    delete _numIntEngine ;
    _numIntEngine= 0;
    if(0 != _numIntegrand) {
      delete _numIntegrand;
      _numIntegrand= 0;
    }
  }

  // All done if there are no arguments to integrate numerically
  if(0 == _intList.getSize()) return kTRUE;
  
  // Bind the appropriate analytic integral (specified by _mode) of our RooRealVar object to
  // those of its arguments that will be integrated out numerically.
  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()) {
    coutE(Integration) << ClassName() << "::" << GetName() << ": failed to create valid integrand." << endl;
    return kFALSE;
  }

  // Create appropriate numeric integrator using factory
  Bool_t isBinned = _function.arg().isBinnedDistribution(_intList) ;
  _numIntEngine = RooNumIntFactory::instance().createIntegrator(*_numIntegrand,*_iconfig,0,isBinned) ;

  if(0 == _numIntEngine || !_numIntEngine->isValid()) {
    coutE(Integration) << ClassName() << "::" << GetName() << ": failed to create valid integrator." << endl;
    return kFALSE;
  }

  cxcoutI(NumIntegration) << "RooRealIntegral::init(" << GetName() << ") using numeric integrator " 
           << _numIntEngine->IsA()->GetName() << " to calculate Int" << _intList << endl ;

  if (_intList.getSize()>3) {
    cxcoutI(NumIntegration) << "RooRealIntegral::init(" << GetName() << ") evaluation requires " << _intList.getSize() << "-D numeric integration step. Evaluation may be slow, sufficient numeric precision for fitting & minimization is not guaranteed" << endl ;
  }

  _restartNumIntEngine = kFALSE ;
  return kTRUE;
}



////////////////////////////////////////////////////////////////////////////////
/// Copy constructor

RooRealIntegral::RooRealIntegral(const RooRealIntegral& other, const char* name) : 
  RooAbsReal(other,name), 
  _valid(other._valid),
  _respectCompSelect(other._respectCompSelect),  
  _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("!func",this,other._function), 
  _iconfig(other._iconfig),
  _sumCat("!sumCat",this,other._sumCat),
  _sumCatIter(0),
  _mode(other._mode),
  _intOperMode(other._intOperMode), 
  _restartNumIntEngine(kFALSE),
  _numIntEngine(0), 
  _numIntegrand(0),
  _rangeName(other._rangeName),
  _params(0),
  _cacheNum(kFALSE)
{
 _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) ;
 }

 other._intList.snapshot(_saveInt) ;
 other._sumList.snapshot(_saveSum) ;

  TRACE_CREATE
}



////////////////////////////////////////////////////////////////////////////////

RooRealIntegral::~RooRealIntegral()
  // Destructor
{
  if (_numIntEngine) delete _numIntEngine ;
  if (_numIntegrand) delete _numIntegrand ;
  if (_funcNormSet) delete _funcNormSet ;
  delete _facListIter ;
  delete _jacListIter ;
  if (_sumCatIter)  delete _sumCatIter ;
  if (_params) delete _params ;

  TRACE_DESTROY
}





////////////////////////////////////////////////////////////////////////////////

RooAbsReal* RooRealIntegral::createIntegral(const RooArgSet& iset, const RooArgSet* nset, const RooNumIntConfig* cfg, const char* rangeName) const 
{
  // Handle special case of no integration with default algorithm
  if (iset.getSize()==0) {
    return RooAbsReal::createIntegral(iset,nset,cfg,rangeName) ;
  }

  // Special handling of integral of integral, return RooRealIntegral that represents integral over all dimensions in one pass
  RooArgSet isetAll(iset) ;
  isetAll.add(_sumList) ;
  isetAll.add(_intList) ;
  isetAll.add(_anaList) ;
  isetAll.add(_facList) ;

  const RooArgSet* newNormSet(0) ;
  RooArgSet* tmp(0) ;
  if (nset && !_funcNormSet) {
    newNormSet = nset ;
  } else if (!nset && _funcNormSet) {
    newNormSet = _funcNormSet ;
  } else if (nset && _funcNormSet) {
    tmp = new RooArgSet ;
    tmp->add(*nset) ;
    tmp->add(*_funcNormSet,kTRUE) ;
    newNormSet = tmp ;
  } 
  RooAbsReal* ret =  _function.arg().createIntegral(isetAll,newNormSet,cfg,rangeName) ;

  if (tmp) {
    delete tmp ;
  }

  return ret ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return value of object. If the cache is clean, return the
/// cached value, otherwise recalculate on the fly and refill
/// the cache

Double_t RooRealIntegral::getValV(const RooArgSet* nset) const
{
//   // fast-track clean-cache processing
//   if (_operMode==AClean) {
//     return _value ;
//   }

  if (nset && nset!=_lastNSet) {
    ((RooAbsReal*) this)->setProxyNormSet(nset) ;    
    _lastNSet = (RooArgSet*) nset ;
  }

  if (isValueOrShapeDirtyAndClear()) {
    _value = traceEval(nset) ;
  } 

  return _value ;
}





////////////////////////////////////////////////////////////////////////////////
/// Perform the integration and return the result

Double_t RooRealIntegral::evaluate() const 
{  

  bool tmp = RooAbsReal::_globalSelectComp;
  if(!_respectCompSelect){
    RooAbsReal::_globalSelectComp = true ;
  }
  
  Double_t retVal(0) ;
  switch (_intOperMode) {    
    
  case Hybrid: 
    {      
      // Cache numeric integrals in >1d expensive object cache
      RooDouble* cacheVal(0) ;
      if ((_cacheNum && _intList.getSize()>0) || _intList.getSize()>=_cacheAllNDim) {
        cacheVal = (RooDouble*) expensiveObjectCache().retrieveObject(GetName(),RooDouble::Class(),parameters())  ;
      }

      if (cacheVal) {
        retVal = *cacheVal ;
   // cout << "using cached value of integral" << GetName() << endl ;
      } else {


        // Find any function dependents that are AClean 
        // and switch them temporarily to ADirty
        Bool_t origState = inhibitDirty() ;
        setDirtyInhibit(kTRUE) ;
        
        // try to initialize our numerical integration engine
        if(!(_valid= initNumIntegrator())) {
          coutE(Integration) << ClassName() << "::" << GetName()
                             << ":evaluate: cannot initialize numerical integrator" << endl;
          RooAbsReal::_globalSelectComp = tmp ;
          return 0;
        }
        
        // Save current integral dependent values 
        _saveInt = _intList ;
        _saveSum = _sumList ;
        
        // Evaluate sum/integral
        retVal = sum() ;
        
        
        // This must happen BEFORE restoring dependents, otherwise no dirty state propagation in restore step
        setDirtyInhibit(origState) ;
        
        // Restore integral dependent values
        _intList=_saveInt ;
        _sumList=_saveSum ;
        
        // Cache numeric integrals in >1d expensive object cache
        if ((_cacheNum && _intList.getSize()>0) || _intList.getSize()>=_cacheAllNDim) {
          RooDouble* val = new RooDouble(retVal) ;
          expensiveObjectCache().registerObject(_function.arg().GetName(),GetName(),*val,parameters())  ;
          //     cout << "### caching value of integral" << GetName() << " in " << &expensiveObjectCache() << endl ;
        }
        
      }
      break ;
    }
  case Analytic:
    {
      retVal =  ((RooAbsReal&)_function.arg()).analyticalIntegralWN(_mode,_funcNormSet,RooNameReg::str(_rangeName)) / jacobianProduct() ;
      cxcoutD(Tracing) << "RooRealIntegral::evaluate_analytic(" << GetName() 
             << ")func = " << _function.arg().IsA()->GetName() << "::" << _function.arg().GetName()
             << " raw = " << retVal << " _funcNormSet = " << (_funcNormSet?*_funcNormSet:RooArgSet()) << endl ;

      
      break ;
    }

  case PassThrough:
    {
      //setDirtyInhibit(kTRUE) ;
      retVal= _function.arg().getVal(_funcNormSet) ;      
      //setDirtyInhibit(kFALSE) ;
      break ;
    }
  }
  

  // Multiply answer with integration ranges of factorized variables
  if (_facList.getSize()>0) {
    RooAbsArg *arg ;
    _facListIter->Reset() ;
    while((arg=(RooAbsArg*)_facListIter->Next())) {
      // Multiply by fit range for 'real' dependents
      if (arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
   RooAbsRealLValue* argLV = (RooAbsRealLValue*)arg ;
   retVal *= (argLV->getMax() - argLV->getMin()) ;
      }
      // Multiply by number of states for category dependents
      if (arg->IsA()->InheritsFrom(RooAbsCategoryLValue::Class())) {
   RooAbsCategoryLValue* argLV = (RooAbsCategoryLValue*)arg ;
   retVal *= argLV->numTypes() ;
      }    
    } 
  }


  if (dologD(Tracing)) {
    cxcoutD(Tracing) << "RooRealIntegral::evaluate(" << GetName() << ") anaInt = " << _anaList << " numInt = " << _intList << _sumList << " mode = " ;
    switch(_mode) {
    case Hybrid: ccoutD(Tracing) << "Hybrid" ; break ;
    case Analytic: ccoutD(Tracing) << "Analytic" ; break ;
    case PassThrough: ccoutD(Tracing) << "PassThrough" ; break ;
    }

    ccxcoutD(Tracing) << "raw*fact = " << retVal << endl ;
  }
  
  RooAbsReal::_globalSelectComp = tmp ;

  return retVal ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return product of jacobian terms originating from analytical integration

Double_t RooRealIntegral::jacobianProduct() const 
{
  if (_jacList.getSize()==0) {
    return 1 ;
  }

  Double_t jacProd(1) ;
  _jacListIter->Reset() ;
  RooAbsRealLValue* arg ;
  while ((arg=(RooAbsRealLValue*)_jacListIter->Next())) {
    jacProd *= arg->jacobian() ;
  }

  // Take fabs() here: if jacobian is negative, min and max are swapped and analytical integral
  // will be positive, so must multiply with positive jacobian.
  return fabs(jacProd) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Perform summation of list of category dependents to be integrated

Double_t RooRealIntegral::sum() const
{  
  if (_sumList.getSize()!=0) {
    // Add integrals for all permutations of categories summed over
    Double_t total(0) ;

    _sumCatIter->Reset() ;
    RooCatType* type ;
    RooSuperCategory* sumCat = (RooSuperCategory*) _sumCat.first() ;
    while((type=(RooCatType*)_sumCatIter->Next())) {
      sumCat->setIndex(type->getVal()) ;
      if (!_rangeName || sumCat->inRange(RooNameReg::str(_rangeName))) {
   total += integrate() / jacobianProduct() ;
      }
    }

    return total ;

  } else {
    // Simply return integral 
    Double_t ret = integrate() / jacobianProduct() ;
    return ret ;
  }
}


////////////////////////////////////////////////////////////////////////////////
/// Perform hybrid numerical/analytical integration over all real-valued dependents

Double_t RooRealIntegral::integrate() const
{
  if (!_numIntEngine) {
    // Trivial case, fully analytical integration
    return ((RooAbsReal&)_function.arg()).analyticalIntegralWN(_mode,_funcNormSet,RooNameReg::str(_rangeName)) ;
  } else {
    return _numIntEngine->calculate()  ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Intercept server redirects and reconfigure internal object accordingly

Bool_t RooRealIntegral::redirectServersHook(const RooAbsCollection& /*newServerList*/, 
                   Bool_t /*mustReplaceAll*/, Bool_t /*nameChange*/, Bool_t /*isRecursive*/) 
{
  _restartNumIntEngine = kTRUE ;

  autoSelectDirtyMode() ;

  // Update contents value caches for _intList and _sumList
  _saveInt.removeAll() ;
  _saveSum.removeAll() ;
  _intList.snapshot(_saveInt) ;
  _sumList.snapshot(_saveSum) ;

  // Delete parameters cache if we have one
  if (_params) {
    delete _params ;
    _params = 0 ;
  }

  return kFALSE ;
}



////////////////////////////////////////////////////////////////////////////////

const RooArgSet& RooRealIntegral::parameters() const
{
  if (!_params) {
    _params = new RooArgSet("params") ;
    
    TIterator* siter = serverIterator() ;
    RooArgSet params ;
    RooAbsArg* server ;
    while((server = (RooAbsArg*)siter->Next())) {
      if (server->isValueServer(*this)) _params->add(*server) ;
    }
    delete siter ;
  }

  return *_params ;
}



////////////////////////////////////////////////////////////////////////////////
/// Dummy

void RooRealIntegral::operModeHook()
{
  if (_operMode==ADirty) {    
//     cout << "RooRealIntegral::operModeHook(" << GetName() << " warning: mode set to ADirty" << endl ;
//     if (TString(GetName()).Contains("FULL")) {
//       cout << "blah" << endl ;
//     }
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Check if current value is valid

Bool_t RooRealIntegral::isValidReal(Double_t /*value*/, Bool_t /*printError*/) const 
{
  return kTRUE ;
}

////////////////////////////////////////////////////////////////////////////////
/// Check if component selection is allowed

Bool_t RooRealIntegral::getAllowComponentSelection() const {
  return _respectCompSelect;
}

////////////////////////////////////////////////////////////////////////////////
/// Set component selection to be allowed/forbidden

void RooRealIntegral::setAllowComponentSelection(Bool_t allow){
  _respectCompSelect = allow;
}

////////////////////////////////////////////////////////////////////////////////
/// Customized printing of arguments of a RooRealIntegral to more intuitively reflect the contents of the
/// integration operation

void RooRealIntegral::printMetaArgs(ostream& os) const
{

  if (intVars().getSize()!=0) {
    os << "Int " ;
  }
  os << _function.arg().GetName() ;
  if (_funcNormSet) {
    os << "_Norm" ;
    os << *_funcNormSet ;
    os << " " ;
  }
 
  // List internally integrated observables and factorizing observables as analytically integrated
  RooArgSet tmp(_anaList) ;
  tmp.add(_facList) ;
  if (tmp.getSize()>0) {
    os << "d[Ana]" ;
    os << tmp ;
    os << " " ;
  }
  
  // List numerically integrated and summed observables as numerically integrated
  RooArgSet tmp2(_intList) ;
  tmp2.add(_sumList) ;
  if (tmp2.getSize()>0) {
    os << " d[Num]" ;
    os << tmp2 ;
    os << " " ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Print the state of this object to the specified output stream.

void RooRealIntegral::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
{
  RooAbsReal::printMultiline(os,contents,verbose,indent) ;
  os << indent << "--- RooRealIntegral ---" << endl; 
  os << indent << "  Integrates ";
  _function.arg().printStream(os,kName|kArgs,kSingleLine,indent);
  TString deeper(indent);
  deeper.Append("  ");
  os << indent << "  operating mode is " 
     << (_intOperMode==Hybrid?"Hybrid":(_intOperMode==Analytic?"Analytic":"PassThrough")) << endl ;
  os << indent << "  Summed discrete args are " << _sumList << endl ;
  os << indent << "  Numerically integrated args are " << _intList << endl;
  os << indent << "  Analytically integrated args using mode " << _mode << " are " << _anaList << endl ;
  os << indent << "  Arguments included in Jacobian are " << _jacList << endl ;
  os << indent << "  Factorized arguments are " << _facList << endl ;
  os << indent << "  Function normalization set " ;
  if (_funcNormSet) 
    _funcNormSet->Print("1") ; 
  else
    os << "<none>" ;
  
  os << endl ;
} 



////////////////////////////////////////////////////////////////////////////////
/// Global switch to cache all integral values that integrate at least ndim dimensions numerically 

void RooRealIntegral::setCacheAllNumeric(Int_t ndim) {
  _cacheAllNDim = ndim ;
}


////////////////////////////////////////////////////////////////////////////////
/// Return minimum dimensions of numeric integration for which values are cached. 

Int_t RooRealIntegral::getCacheAllNumeric() 
{
  return _cacheAllNDim ;
}