RooAbsReal.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)             *
 *****************************************************************************/

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

/** \class RooAbsReal

   RooAbsReal is the common abstract base class for objects that represent a
   real value and implements functionality common to all real-valued objects
   such as the ability to plot them, to construct integrals of them, the
   ability to advertise (partial) analytical integrals etc.

   Implementation of RooAbsReal may be derived, thus no interface
   is provided to modify the contents.

   \ingroup Roofitcore
*/




#include <sys/types.h>


#include "RooFit.h"
#include "RooMsgService.h"

#include "RooAbsReal.h"
#include "RooAbsReal.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooBinning.h"
#include "RooPlot.h"
#include "RooCurve.h"
#include "RooRealVar.h"
#include "RooArgProxy.h"
#include "RooFormulaVar.h"
#include "RooRealBinding.h"
#include "RooRealIntegral.h"
#include "RooAbsCategoryLValue.h"
#include "RooCustomizer.h"
#include "RooAbsData.h"
#include "RooScaledFunc.h"
#include "RooAddPdf.h"
#include "RooCmdConfig.h"
#include "RooCategory.h"
#include "RooNumIntConfig.h"
#include "RooAddition.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooDataWeightedAverage.h"
#include "RooNumRunningInt.h"
#include "RooGlobalFunc.h"
#include "RooParamBinning.h"
#include "RooProfileLL.h"
#include "RooFunctor.h"
#include "RooDerivative.h"
#include "RooGenFunction.h"
#include "RooMultiGenFunction.h"
#include "RooCmdConfig.h"
#include "RooXYChi2Var.h"
#include "RooMinuit.h"
#ifndef __ROOFIT_NOROOMINIMIZER
#include "RooMinimizer.h"
#endif
#include "RooChi2Var.h"
#include "RooFitResult.h"
#include "RooAbsMoment.h"
#include "RooMoment.h"
#include "RooFirstMoment.h"
#include "RooSecondMoment.h"
#include "RooBrentRootFinder.h"
#include "RooVectorDataStore.h"
#include "RooCachedReal.h"

#include "Riostream.h"

#include "Math/IFunction.h"
#include "TMath.h"
#include "TObjString.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TAttLine.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TMatrixD.h"
#include "TVector.h"

#include <sstream>

using namespace std ;

ClassImp(RooAbsReal)
;

Bool_t RooAbsReal::_cacheCheck(kFALSE) ;
Bool_t RooAbsReal::_globalSelectComp = kFALSE ;
Bool_t RooAbsReal::_hideOffset = kTRUE ;

void RooAbsReal::setHideOffset(Bool_t flag) { _hideOffset = flag ; }
Bool_t RooAbsReal::hideOffset() { return _hideOffset ; }

RooAbsReal::ErrorLoggingMode RooAbsReal::_evalErrorMode = RooAbsReal::PrintErrors ;
Int_t RooAbsReal::_evalErrorCount = 0 ;
map<const RooAbsArg*,pair<string,list<RooAbsReal::EvalError> > > RooAbsReal::_evalErrorList ;


////////////////////////////////////////////////////////////////////////////////
/// coverity[UNINIT_CTOR]
/// Default constructor

RooAbsReal::RooAbsReal() : _specIntegratorConfig(0), _treeVar(kFALSE), _selectComp(kTRUE), _lastNSet(0)
{
}



////////////////////////////////////////////////////////////////////////////////
/// Constructor with unit label

RooAbsReal::RooAbsReal(const char *name, const char *title, const char *unit) :
  RooAbsArg(name,title), _plotMin(0), _plotMax(0), _plotBins(100),
  _value(0),  _unit(unit), _forceNumInt(kFALSE), _specIntegratorConfig(0), _treeVar(kFALSE), _selectComp(kTRUE), _lastNSet(0)
{
  setValueDirty() ;
  setShapeDirty() ;

}



////////////////////////////////////////////////////////////////////////////////
/// Constructor with plot range and unit label

RooAbsReal::RooAbsReal(const char *name, const char *title, Double_t inMinVal,
             Double_t inMaxVal, const char *unit) :
  RooAbsArg(name,title), _plotMin(inMinVal), _plotMax(inMaxVal), _plotBins(100),
  _value(0), _unit(unit), _forceNumInt(kFALSE), _specIntegratorConfig(0), _treeVar(kFALSE), _selectComp(kTRUE), _lastNSet(0)
{
  setValueDirty() ;
  setShapeDirty() ;

}



////////////////////////////////////////////////////////////////////////////////
/// coverity[UNINIT_CTOR]
/// Copy constructor

RooAbsReal::RooAbsReal(const RooAbsReal& other, const char* name) :
  RooAbsArg(other,name), _plotMin(other._plotMin), _plotMax(other._plotMax),
  _plotBins(other._plotBins), _value(other._value), _unit(other._unit), _label(other._label),
  _forceNumInt(other._forceNumInt), _treeVar(other._treeVar), _selectComp(other._selectComp), _lastNSet(0)
{
  if (other._specIntegratorConfig) {
    _specIntegratorConfig = new RooNumIntConfig(*other._specIntegratorConfig) ;
  } else {
    _specIntegratorConfig = 0 ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Destructor

RooAbsReal::~RooAbsReal()
{
  if (_specIntegratorConfig) delete _specIntegratorConfig ;
}



////////////////////////////////////////////////////////////////////////////////
/// Equality operator comparing to a Double_t

Bool_t RooAbsReal::operator==(Double_t value) const
{
  return (getVal()==value) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Equality operator when comparing to another RooAbsArg.
/// Only functional when the other arg is a RooAbsReal

Bool_t RooAbsReal::operator==(const RooAbsArg& other)
{
  const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
  return otherReal ? operator==(otherReal->getVal()) : kFALSE ;
}


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

Bool_t RooAbsReal::isIdentical(const RooAbsArg& other, Bool_t assumeSameType)
{
  if (!assumeSameType) {
    const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
    return otherReal ? operator==(otherReal->getVal()) : kFALSE ;
  } else {
    return getVal()==((RooAbsReal&)other).getVal() ;
  }
}


////////////////////////////////////////////////////////////////////////////////
/// Return this variable's title string. If appendUnit is true and
/// this variable has units, also append a string " (<unit>)".

TString RooAbsReal::getTitle(Bool_t appendUnit) const
{
  TString title(GetTitle());
  if(appendUnit && 0 != strlen(getUnit())) {
    title.Append(" (");
    title.Append(getUnit());
    title.Append(")");
  }
  return title;
}



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

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

  if (isValueDirtyAndClear()) {
    _value = traceEval(nset) ;
    //     clearValueDirty() ;
  }
  //   cout << "RooAbsReal::getValV(" << GetName() << ") writing _value = " << _value << endl ;

  Double_t ret(_value) ;
  if (hideOffset()) ret += offset() ;

  return ret ;
}


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

Int_t RooAbsReal::numEvalErrorItems()
{
  return _evalErrorList.size() ;
}


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

RooAbsReal::EvalErrorIter RooAbsReal::evalErrorIter()
{
  return _evalErrorList.begin() ;
}


////////////////////////////////////////////////////////////////////////////////
/// Calculate current value of object, with error tracing wrapper

Double_t RooAbsReal::traceEval(const RooArgSet* /*nset*/) const
{
  Double_t value = evaluate() ;

  if (TMath::IsNaN(value)) {
    logEvalError("function value is NAN") ;
  }

  //cxcoutD(Tracing) << "RooAbsReal::getValF(" << GetName() << ") operMode = " << _operMode << " recalculated, new value = " << value << endl ;

  //Standard tracing code goes here
  if (!isValidReal(value)) {
    coutW(Tracing) << "RooAbsReal::traceEval(" << GetName()
         << "): validation failed: " << value << endl ;
  }

  //Call optional subclass tracing code
  //   traceEvalHook(value) ;

  return value ;
}



////////////////////////////////////////////////////////////////////////////////
/// Variant of getAnalyticalIntegral that is also passed the normalization set
/// that should be applied to the integrand of which the integral is request.
/// For certain operator p.d.f it is useful to overload this function rather
/// than analyticalIntegralWN() as the additional normalization information
/// may be useful in determining a more efficient decomposition of the
/// requested integral

Int_t RooAbsReal::getAnalyticalIntegralWN(RooArgSet& allDeps, RooArgSet& analDeps,
                 const RooArgSet* /*normSet*/, const char* rangeName) const
{
  return _forceNumInt ? 0 : getAnalyticalIntegral(allDeps,analDeps,rangeName) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Interface function getAnalyticalIntergral advertises the
/// analytical integrals that are supported. 'integSet'
/// is the set of dependents for which integration is requested. The
/// function should copy the subset of dependents it can analytically
/// integrate to anaIntSet and return a unique identification code for
/// this integration configuration.  If no integration can be
/// performed, zero should be returned.

Int_t RooAbsReal::getAnalyticalIntegral(RooArgSet& /*integSet*/, RooArgSet& /*anaIntSet*/, const char* /*rangeName*/) const
{
  return 0 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Implements the actual analytical integral(s) advertised by
/// getAnalyticalIntegral.  This functions will only be called with
/// codes returned by getAnalyticalIntegral, except code zero.

Double_t RooAbsReal::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
{
//   cout << "RooAbsReal::analyticalIntegralWN(" << GetName() << ") code = " << code << " normSet = " << (normSet?*normSet:RooArgSet()) << endl ;
  if (code==0) return getVal(normSet) ;
  return analyticalIntegral(code,rangeName) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Implements the actual analytical integral(s) advertised by
/// getAnalyticalIntegral.  This functions will only be called with
/// codes returned by getAnalyticalIntegral, except code zero.

Double_t RooAbsReal::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
{
  // By default no analytical integrals are implemented
  coutF(Eval)  << "RooAbsReal::analyticalIntegral(" << GetName() << ") code " << code << " not implemented" << endl ;
  return 0 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Get the label associated with the variable

const char *RooAbsReal::getPlotLabel() const
{
  return _label.IsNull() ? fName.Data() : _label.Data();
}



////////////////////////////////////////////////////////////////////////////////
/// Set the label associated with this variable

void RooAbsReal::setPlotLabel(const char *label)
{
  _label= label;
}



////////////////////////////////////////////////////////////////////////////////
///Read object contents from stream (dummy for now)

Bool_t RooAbsReal::readFromStream(istream& /*is*/, Bool_t /*compact*/, Bool_t /*verbose*/)
{
  return kFALSE ;
}



////////////////////////////////////////////////////////////////////////////////
///Write object contents to stream (dummy for now)

void RooAbsReal::writeToStream(ostream& /*os*/, Bool_t /*compact*/) const
{
}



////////////////////////////////////////////////////////////////////////////////
/// Print object value

void RooAbsReal::printValue(ostream& os) const
{
  os << getVal() ;
}



////////////////////////////////////////////////////////////////////////////////
/// Structure printing

void RooAbsReal::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
{
  RooAbsArg::printMultiline(os,contents,verbose,indent) ;
  os << indent << "--- RooAbsReal ---" << endl;
  TString unit(_unit);
  if(!unit.IsNull()) unit.Prepend(' ');
  //os << indent << "  Value = " << getVal() << unit << endl;
  os << endl << indent << "  Plot label is \"" << getPlotLabel() << "\"" << endl;

}


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

Bool_t RooAbsReal::isValid() const
{
  return isValidReal(_value) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Interface function to check if given value is a valid value for this object.
/// This default implementation considers all values valid

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




////////////////////////////////////////////////////////////////////////////////
/// Create a RooProfileLL object that eliminates all nuisance parameters in the
/// present function. The nuisance parameters are defined as all parameters
/// of the function except the stated paramsOfInterest

RooAbsReal* RooAbsReal::createProfile(const RooArgSet& paramsOfInterest)
{
  // Construct name of profile object
  TString name(Form("%s_Profile[",GetName())) ;
  TIterator* iter = paramsOfInterest.createIterator() ;
  RooAbsArg* arg ;
  Bool_t first(kTRUE) ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (first) {
      first=kFALSE ;
    } else {
      name.Append(",") ;
    }
    name.Append(arg->GetName()) ;
  }
  delete iter ;
  name.Append("]") ;

  // Create and return profile object
  return new RooProfileLL(name.Data(),Form("Profile of %s",GetTitle()),*this,paramsOfInterest) ;
}






////////////////////////////////////////////////////////////////////////////////
/// Create an object that represents the integral of the function over one or more observables listed in iset
/// The actual integration calculation is only performed when the return object is evaluated. The name
/// of the integral object is automatically constructed from the name of the input function, the variables
/// it integrates and the range integrates over
///
/// The following named arguments are accepted
///
/// NormSet(const RooArgSet&)            -- Specify normalization set, mostly useful when working with PDFS
/// NumIntConfig(const RooNumIntConfig&) -- Use given configuration for any numeric integration, if necessary
/// Range(const char* name)              -- Integrate only over given range. Multiple ranges may be specified
///                                         by passing multiple Range() arguments

RooAbsReal* RooAbsReal::createIntegral(const RooArgSet& iset, const RooCmdArg& arg1, const RooCmdArg& arg2,
                   const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
                   const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
{


  // Define configuration for this method
  RooCmdConfig pc(Form("RooAbsReal::createIntegral(%s)",GetName())) ;
  pc.defineString("rangeName","RangeWithName",0,"",kTRUE) ;
  pc.defineObject("normSet","NormSet",0,0) ;
  pc.defineObject("numIntConfig","NumIntConfig",0,0) ;

  // Process & check varargs
  pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
  if (!pc.ok(kTRUE)) {
    return 0 ;
  }

  // Extract values from named arguments
  const char* rangeName = pc.getString("rangeName",0,kTRUE) ;
  const RooArgSet* nset = static_cast<const RooArgSet*>(pc.getObject("normSet",0)) ;
  const RooNumIntConfig* cfg = static_cast<const RooNumIntConfig*>(pc.getObject("numIntConfig",0)) ;

  return createIntegral(iset,nset,cfg,rangeName) ;
}





////////////////////////////////////////////////////////////////////////////////
/// Create an object that represents the integral of the function over one or more observables listed in iset
/// The actual integration calculation is only performed when the return object is evaluated. The name
/// of the integral object is automatically constructed from the name of the input function, the variables
/// it integrates and the range integrates over. If nset is specified the integrand is request
/// to be normalized over nset (only meaningful when the integrand is a pdf). If rangename is specified
/// the integral is performed over the named range, otherwise it is performed over the domain of each
/// integrated observable. If cfg is specified it will be used to configure any numeric integration
/// aspect of the integral. It will not force the integral to be performed numerically, which is
/// decided automatically by RooRealIntegral

RooAbsReal* RooAbsReal::createIntegral(const RooArgSet& iset, const RooArgSet* nset,
                   const RooNumIntConfig* cfg, const char* rangeName) const
{
  if (!rangeName || strchr(rangeName,',')==0) {
    // Simple case: integral over full range or single limited range
    return createIntObj(iset,nset,cfg,rangeName) ;
  }

  // Integral over multiple ranges
  RooArgSet components ;

  TObjArray* oa = TString(rangeName).Tokenize(",");

  for( Int_t i=0; i < oa->GetEntries(); ++i) {
    TObjString* os = (TObjString*) (*oa)[i];
    if(!os) break;
    RooAbsReal* compIntegral = createIntObj(iset,nset,cfg,os->GetString().Data()) ;
    components.add(*compIntegral) ;
  }
  delete oa;

  TString title(GetTitle()) ;
  title.Prepend("Integral of ") ;
  TString fullName(GetName()) ;
  fullName.Append(integralNameSuffix(iset,nset,rangeName)) ;

  return new RooAddition(fullName.Data(),title.Data(),components,kTRUE) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for createIntegral that creates the actual integreal object

RooAbsReal* RooAbsReal::createIntObj(const RooArgSet& iset2, const RooArgSet* nset2,
                 const RooNumIntConfig* cfg, const char* rangeName) const
{
  // Make internal use copies of iset and nset
  RooArgSet iset(iset2) ;
  const RooArgSet* nset = nset2 ;


  // Initialize local variables perparing for recursive loop
  Bool_t error = kFALSE ;
  const RooAbsReal* integrand = this ;
  RooAbsReal* integral = 0 ;

  // Handle trivial case of no integration here explicitly
  if (iset.getSize()==0) {

    TString title(GetTitle()) ;
    title.Prepend("Integral of ") ;

    TString name(GetName()) ;
    name.Append(integralNameSuffix(iset,nset,rangeName)) ;

    return new RooRealIntegral(name,title,*this,iset,nset,cfg,rangeName) ;
  }

  // Process integration over remaining integration variables
  while(iset.getSize()>0) {


    // Find largest set of observables that can be integrated in one go
    RooArgSet innerSet ;
    findInnerMostIntegration(iset,innerSet,rangeName) ;

    // If largest set of observables that can be integrated is empty set, problem was ill defined
    // Postpone error messaging and handling to end of function, exit loop here
    if (innerSet.getSize()==0) {
      error = kTRUE ;
      break ;
    }

    // Prepare name and title of integral to be created
    TString title(integrand->GetTitle()) ;
    title.Prepend("Integral of ") ;

    TString name(integrand->GetName()) ;
    name.Append(integrand->integralNameSuffix(innerSet,nset,rangeName)) ;

    // Construct innermost integral
    integral = new RooRealIntegral(name,title,*integrand,innerSet,nset,cfg,rangeName) ;

    // Integral of integral takes ownership of innermost integral
    if (integrand != this) {
      integral->addOwnedComponents(*integrand) ;
    }

    // Remove already integrated observables from to-do list
    iset.remove(innerSet) ;

    // Send info message on recursion if needed
    if (integrand == this && iset.getSize()>0) {
      coutI(Integration) << GetName() << " : multidimensional integration over observables with parameterized ranges in terms of other integrated observables detected, using recursive integration strategy to construct final integral" << endl ;
    }

    // Prepare for recursion, next integral should integrate last integrand
    integrand = integral ;


    // Only need normalization set in innermost integration
    nset = 0 ;
  }

  if (error) {
    coutE(Integration) << GetName() << " : ERROR while defining recursive integral over observables with parameterized integration ranges, please check that integration rangs specify uniquely defined integral " << endl;
    delete integral ;
    integral = 0 ;
    return integral ;
  }


  // After-burner: apply interpolating cache on (numeric) integral if requested by user
  const char* cacheParamsStr = getStringAttribute("CACHEPARAMINT") ;
  if (cacheParamsStr && strlen(cacheParamsStr)) {

    RooArgSet* intParams = integral->getVariables() ;

    RooNameSet cacheParamNames ;
    cacheParamNames.setNameList(cacheParamsStr) ;
    RooArgSet* cacheParams = cacheParamNames.select(*intParams) ;

    if (cacheParams->getSize()>0) {
      cxcoutD(Caching) << "RooAbsReal::createIntObj(" << GetName() << ") INFO: constructing " << cacheParams->getSize()
           << "-dim value cache for integral over " << iset2 << " as a function of " << *cacheParams << " in range " << (rangeName?rangeName:"<none>") <<  endl ;
      string name = Form("%s_CACHE_[%s]",integral->GetName(),cacheParams->contentsString().c_str()) ;
      RooCachedReal* cachedIntegral = new RooCachedReal(name.c_str(),name.c_str(),*integral,*cacheParams) ;
      cachedIntegral->setInterpolationOrder(2) ;
      cachedIntegral->addOwnedComponents(*integral) ;
      cachedIntegral->setCacheSource(kTRUE) ;
      if (integral->operMode()==ADirty) {
   cachedIntegral->setOperMode(ADirty) ;
      }
      //cachedIntegral->disableCache(kTRUE) ;
      integral = cachedIntegral ;
    }

    delete cacheParams ;
    delete intParams ;
  }

  return integral ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for createIntObj() that aids in the construct of recursive integrals
/// over functions with multiple observables with parameterized ranges. This function
/// finds in a given set allObs over which integration is requested the largeset subset
/// of observables that can be integrated simultaneously. This subset consists of
/// observables with fixed ranges and observables with parameterized ranges whose
/// parameterization does not depend on any observable that is also integrated.

void RooAbsReal::findInnerMostIntegration(const RooArgSet& allObs, RooArgSet& innerObs, const char* rangeName) const
{
  // Make lists of
  // a) integrated observables with fixed ranges,
  // b) integrated observables with parameterized ranges depending on other integrated observables
  // c) integrated observables used in definition of any parameterized ranges of integrated observables
  RooArgSet obsWithFixedRange(allObs) ;
  RooArgSet obsWithParamRange ;
  RooArgSet obsServingAsRangeParams ;

  // Loop over all integrated observables
  TIterator* oiter = allObs.createIterator() ;
  RooAbsArg* aarg ;
  while((aarg=(RooAbsArg*)oiter->Next())) {
    // Check if observable is real-valued lvalue
    RooAbsRealLValue* arglv = dynamic_cast<RooAbsRealLValue*>(aarg) ;
    if (arglv) {

      // Check if range is parameterized
      RooAbsBinning& binning = arglv->getBinning(rangeName,kFALSE,kTRUE) ;
      if (binning.isParameterized()) {
   RooArgSet* loBoundObs = binning.lowBoundFunc()->getObservables(allObs) ;
   RooArgSet* hiBoundObs = binning.highBoundFunc()->getObservables(allObs) ;

   // Check if range parameterization depends on other integrated observables
   if (loBoundObs->overlaps(allObs) || hiBoundObs->overlaps(allObs)) {
     obsWithParamRange.add(*aarg) ;
     obsWithFixedRange.remove(*aarg) ;
     obsServingAsRangeParams.add(*loBoundObs,kFALSE) ;
     obsServingAsRangeParams.add(*hiBoundObs,kFALSE) ;
   }
   delete loBoundObs ;
   delete hiBoundObs ;
      }
    }
  }
  delete oiter ;

  // Make list of fixed-range observables that are _not_ involved in the parameterization of ranges of other observables
  RooArgSet obsWithFixedRangeNP(obsWithFixedRange) ;
  obsWithFixedRangeNP.remove(obsServingAsRangeParams) ;

  // Make list of param-range observables that are _not_ involved in the parameterization of ranges of other observables
  RooArgSet obsWithParamRangeNP(obsWithParamRange) ;
  obsWithParamRangeNP.remove(obsServingAsRangeParams) ;

  // Construct inner-most integration: over observables (with fixed or param range) not used in any other param range definitions
  innerObs.removeAll() ;
  innerObs.add(obsWithFixedRangeNP) ;
  innerObs.add(obsWithParamRangeNP) ;

}


////////////////////////////////////////////////////////////////////////////////
/// Construct string with unique suffix name to give to integral object that encodes
/// integrated observables, normalization observables and the integration range name

TString RooAbsReal::integralNameSuffix(const RooArgSet& iset, const RooArgSet* nset, const char* rangeName, Bool_t omitEmpty) const
{
  TString name ;
  if (iset.getSize()>0) {

    RooArgSet isetTmp(iset) ;
    isetTmp.sort() ;

    name.Append("_Int[") ;
    TIterator* iter = isetTmp.createIterator() ;
    RooAbsArg* arg ;
    Bool_t first(kTRUE) ;
    while((arg=(RooAbsArg*)iter->Next())) {
      if (first) {
   first=kFALSE ;
      } else {
   name.Append(",") ;
      }
      name.Append(arg->GetName()) ;
    }
    delete iter ;
    if (rangeName) {
      name.Append("|") ;
      name.Append(rangeName) ;
    }
    name.Append("]");
  } else if (!omitEmpty) {
    name.Append("_Int[]") ;
  }

  if (nset && nset->getSize()>0 ) {

    RooArgSet nsetTmp(*nset) ;
    nsetTmp.sort() ;

    name.Append("_Norm[") ;
    Bool_t first(kTRUE);
    TIterator* iter  = nsetTmp.createIterator() ;
    RooAbsArg* arg ;
    while((arg=(RooAbsArg*)iter->Next())) {
      if (first) {
   first=kFALSE ;
      } else {
   name.Append(",") ;
      }
      name.Append(arg->GetName()) ;
    }
    delete iter ;
    const RooAbsPdf* thisPdf = dynamic_cast<const RooAbsPdf*>(this) ;
    if (thisPdf && thisPdf->normRange()) {
      name.Append("|") ;
      name.Append(thisPdf->normRange()) ;
    }
    name.Append("]") ;
  }

  return name ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for plotOn() that creates a projection of a function or p.d.f
/// to be plotted on a RooPlot.

const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars,
                                               RooArgSet*& cloneSet) const
{
  return createPlotProjection(depVars,&projVars,cloneSet) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for plotOn() that creates a projection of a function or p.d.f
/// to be plotted on a RooPlot.

const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars) const
{
  RooArgSet* cloneSet = new RooArgSet() ;
  return createPlotProjection(depVars,&projVars,cloneSet) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for plotOn() that creates a projection of a function or p.d.f
/// to be plotted on a RooPlot.
///
/// Create a new object G that represents the normalized projection:
///
///             Integral [ F[x,y,p] , { y } ]
///  G[x,p] = ---------------------------------
///            Integral [ F[x,y,p] , { x,y } ]
///
/// where F[x,y,p] is the function we represent, "x" are the
/// specified dependentVars, "y" are the specified projectedVars, and
/// "p" are our remaining variables ("parameters"). Return a
/// pointer to the newly created object, or else zero in case of an
/// error.  The caller is responsible for deleting the contents of
/// cloneSet (which includes the returned projection object)

const RooAbsReal *RooAbsReal::createPlotProjection(const RooArgSet &dependentVars, const RooArgSet *projectedVars,
                     RooArgSet *&cloneSet, const char* rangeName, const RooArgSet* condObs) const
{
  // Get the set of our leaf nodes
  RooArgSet leafNodes;
  RooArgSet treeNodes;
  leafNodeServerList(&leafNodes,this);
  treeNodeServerList(&treeNodes,this) ;


  // Check that the dependents are all fundamental. Filter out any that we
  // do not depend on, and make substitutions by name in our leaf list.
  // Check for overlaps with the projection variables.

  TIterator *dependentIterator= dependentVars.createIterator();
  assert(0 != dependentIterator);
  const RooAbsArg *arg = 0;
  while((arg= (const RooAbsArg*)dependentIterator->Next())) {
    if(!arg->isFundamental() && !dynamic_cast<const RooAbsLValue*>(arg)) {
      coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: variable \"" << arg->GetName()
      << "\" of wrong type: " << arg->ClassName() << endl;
      delete dependentIterator;
      return 0;
    }

    RooAbsArg *found= treeNodes.find(arg->GetName());
    if(!found) {
      coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
            << "\" is not a dependent and will be ignored." << endl;
      continue;
    }
    if(found != arg) {
      if (leafNodes.find(found->GetName())) {
   leafNodes.replace(*found,*arg);
      } else {
   leafNodes.add(*arg) ;

   // Remove any dependents of found, replace by dependents of LV node
   RooArgSet* lvDep = arg->getObservables(&leafNodes) ;
   RooAbsArg* lvs ;
   TIterator* iter = lvDep->createIterator() ;
   while((lvs=(RooAbsArg*)iter->Next())) {
     RooAbsArg* tmp = leafNodes.find(lvs->GetName()) ;
     if (tmp) {
       leafNodes.remove(*tmp) ;
       leafNodes.add(*lvs) ;
     }
   }
   delete iter ;

      }
    }

    // check if this arg is also in the projection set
    if(0 != projectedVars && projectedVars->find(arg->GetName())) {
      coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
            << "\" cannot be both a dependent and a projected variable." << endl;
      delete dependentIterator;
      return 0;
    }
  }

  // Remove the projected variables from the list of leaf nodes, if necessary.
  if(0 != projectedVars) leafNodes.remove(*projectedVars,kTRUE);

  // Make a deep-clone of ourself so later operations do not disturb our original state
  cloneSet= (RooArgSet*)RooArgSet(*this).snapshot(kTRUE);
  if (!cloneSet) {
    coutE(Plotting) << "RooAbsPdf::createPlotProjection(" << GetName() << ") Couldn't deep-clone PDF, abort," << endl ;
    return 0 ;
  }
  RooAbsReal *theClone= (RooAbsReal*)cloneSet->find(GetName());

  // The remaining entries in our list of leaf nodes are the the external
  // dependents (x) and parameters (p) of the projection. Patch them back
  // into the theClone. This orphans the nodes they replace, but the orphans
  // are still in the cloneList and so will be cleaned up eventually.
  //cout << "redirection leafNodes : " ; leafNodes.Print("1") ;

  RooArgSet* plotLeafNodes = (RooArgSet*) leafNodes.selectCommon(dependentVars) ;
  theClone->recursiveRedirectServers(*plotLeafNodes,kFALSE,kFALSE,kFALSE);
  delete plotLeafNodes ;

  // Create the set of normalization variables to use in the projection integrand
  RooArgSet normSet(dependentVars);
  if(0 != projectedVars) normSet.add(*projectedVars);
  if(0 != condObs) {
    normSet.remove(*condObs,kTRUE,kTRUE) ;
  }

  // Try to create a valid projection integral. If no variables are to be projected,
  // create a null projection anyway to bind our normalization over the dependents
  // consistently with the way they would be bound with a non-trivial projection.
  RooArgSet empty;
  if(0 == projectedVars) projectedVars= &empty;

  TString name = GetName() ;
  name += integralNameSuffix(*projectedVars,&normSet,rangeName,kTRUE) ;

  TString title(GetTitle());
  title.Prepend("Projection of ");


  RooAbsReal* projected= theClone->createIntegral(*projectedVars,normSet,rangeName) ;

  if(0 == projected || !projected->isValid()) {
    coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: cannot integrate out ";
    projectedVars->printStream(cout,kName|kArgs,kSingleLine);
    // cleanup and exit
    if(0 != projected) delete projected;
    delete dependentIterator;
    return 0;
  }

  projected->SetName(name.Data()) ;
  projected->SetTitle(title.Data()) ;

  // Add the projection integral to the cloneSet so that it eventually gets cleaned up by the caller.
  cloneSet->addOwned(*projected);

  // cleanup
  delete dependentIterator;

  // return a const pointer to remind the caller that they do not delete the returned object
  // directly (it is contained in the cloneSet instead).
  return projected;
}




////////////////////////////////////////////////////////////////////////////////
/// Fill the ROOT histogram 'hist' with values sampled from this
/// function at the bin centers.  Our value is calculated by first
/// integrating out any variables in projectedVars and then scaling
/// the result by scaleFactor. Returns a pointer to the input
/// histogram, or zero in case of an error. The input histogram can
/// be any TH1 subclass, and therefore of arbitrary
/// dimension. Variables are matched with the (x,y,...) dimensions of
/// the input histogram according to the order in which they appear
/// in the input plotVars list. If scaleForDensity is true the
/// histogram is filled with a the functions density rather than
/// the functions value (i.e. the value at the bin center is multiplied
/// with bin volume)

TH1 *RooAbsReal::fillHistogram(TH1 *hist, const RooArgList &plotVars,
                Double_t scaleFactor, const RooArgSet *projectedVars, Bool_t scaleForDensity,
                const RooArgSet* condObs, Bool_t setError) const
{
  // Do we have a valid histogram to use?
  if(0 == hist) {
    coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: no valid histogram to fill" << endl;
    return 0;
  }

  // Check that the number of plotVars matches the input histogram's dimension
  Int_t hdim= hist->GetDimension();
  if(hdim != plotVars.getSize()) {
    coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: plotVars has the wrong dimension" << endl;
    return 0;
  }


  // Check that the plot variables are all actually RooRealVars and print a warning if we do not
  // explicitly depend on one of them. Fill a set (not list!) of cloned plot variables.
  RooArgSet plotClones;
  for(Int_t index= 0; index < plotVars.getSize(); index++) {
    const RooAbsArg *var= plotVars.at(index);
    const RooRealVar *realVar= dynamic_cast<const RooRealVar*>(var);
    if(0 == realVar) {
      coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot plot variable \"" << var->GetName()
      << "\" of type " << var->ClassName() << endl;
      return 0;
    }
    if(!this->dependsOn(*realVar)) {
      coutE(InputArguments) << ClassName() << "::" << GetName()
      << ":fillHistogram: WARNING: variable is not an explicit dependent: " << realVar->GetName() << endl;
    }
    plotClones.addClone(*realVar,kTRUE); // do not complain about duplicates
  }

  // Reconnect all plotClones to each other, imported when plotting N-dim integrals with entangled parameterized ranges
  TIterator* pciter= plotClones.createIterator() ;
  RooAbsArg* pc ;
  while((pc=(RooAbsArg*)pciter->Next())) {
    pc->recursiveRedirectServers(plotClones,kFALSE,kFALSE,kTRUE) ;
  }

  delete pciter ;

  // Call checkObservables
  RooArgSet allDeps(plotClones) ;
  if (projectedVars) {
    allDeps.add(*projectedVars) ;
  }
  if (checkObservables(&allDeps)) {
    coutE(InputArguments) << "RooAbsReal::fillHistogram(" << GetName() << ") error in checkObservables, abort" << endl ;
    return hist ;
  }

  // Create a standalone projection object to use for calculating bin contents
  RooArgSet *cloneSet = 0;
  const RooAbsReal *projected= createPlotProjection(plotClones,projectedVars,cloneSet,0,condObs);

  cxcoutD(Plotting) << "RooAbsReal::fillHistogram(" << GetName() << ") plot projection object is " << projected->GetName() << endl ;

  // Prepare to loop over the histogram bins
  Int_t xbins(0),ybins(1),zbins(1);
  RooRealVar *xvar = 0;
  RooRealVar *yvar = 0;
  RooRealVar *zvar = 0;
  TAxis *xaxis = 0;
  TAxis *yaxis = 0;
  TAxis *zaxis = 0;
  switch(hdim) {
  case 3:
    zbins= hist->GetNbinsZ();
    zvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(2)->GetName()));
    zaxis= hist->GetZaxis();
    assert(0 != zvar && 0 != zaxis);
    if (scaleForDensity) {
      scaleFactor*= (zaxis->GetXmax() - zaxis->GetXmin())/zbins;
    }
    // fall through to next case...
  case 2:
    ybins= hist->GetNbinsY();
    yvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(1)->GetName()));
    yaxis= hist->GetYaxis();
    assert(0 != yvar && 0 != yaxis);
    if (scaleForDensity) {
      scaleFactor*= (yaxis->GetXmax() - yaxis->GetXmin())/ybins;
    }
    // fall through to next case...
  case 1:
    xbins= hist->GetNbinsX();
    xvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(0)->GetName()));
    xaxis= hist->GetXaxis();
    assert(0 != xvar && 0 != xaxis);
    if (scaleForDensity) {
      scaleFactor*= (xaxis->GetXmax() - xaxis->GetXmin())/xbins;
    }
    break;
  default:
    coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot fill histogram with "
           << hdim << " dimensions" << endl;
    break;
  }

  // Loop over the input histogram's bins and fill each one with our projection's
  // value, calculated at the center.
  RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
  Int_t xbin(0),ybin(0),zbin(0);
  Int_t bins= xbins*ybins*zbins;
  for(Int_t bin= 0; bin < bins; bin++) {
    switch(hdim) {
    case 3:
      if(bin % (xbins*ybins) == 0) {
   zbin++;
   zvar->setVal(zaxis->GetBinCenter(zbin));
      }
      // fall through to next case...
    case 2:
      if(bin % xbins == 0) {
   ybin= (ybin%ybins) + 1;
   yvar->setVal(yaxis->GetBinCenter(ybin));
      }
      // fall through to next case...
    case 1:
      xbin= (xbin%xbins) + 1;
      xvar->setVal(xaxis->GetBinCenter(xbin));
      break;
    default:
      coutE(InputArguments) << "RooAbsReal::fillHistogram: Internal Error!" << endl;
      break;
    }

    Double_t result= scaleFactor*projected->getVal();
    if (RooAbsReal::numEvalErrors()>0) {
      coutW(Plotting) << "WARNING: Function evaluation error(s) at coordinates [x]=" << xvar->getVal() ;
      if (hdim==2) ccoutW(Plotting) << " [y]=" << yvar->getVal() ;
      if (hdim==3) ccoutW(Plotting) << " [z]=" << zvar->getVal() ;
      ccoutW(Plotting) << endl ;
      // RooAbsReal::printEvalErrors(ccoutW(Plotting),10) ;
      result = 0 ;
    }
    RooAbsReal::clearEvalErrorLog() ;

    hist->SetBinContent(hist->GetBin(xbin,ybin,zbin),result);
    if (setError) {
      hist->SetBinError(hist->GetBin(xbin,ybin,zbin),sqrt(result)) ;
    }

    //cout << "bin " << bin << " -> (" << xbin << "," << ybin << "," << zbin << ") = " << result << endl;
  }
  RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;

  // cleanup
  delete cloneSet;

  return hist;
}



////////////////////////////////////////////////////////////////////////////////
/// Fill a RooDataHist with values sampled from this function at the
/// bin centers.  If extendedMode is true, the p.d.f. values is multiplied
/// by the number of expected events in each bin
///
/// An optional scaling by a given scaleFactor can be performed.
/// Returns a pointer to the input RooDataHist, or zero
/// in case of an error.
///
/// If correctForBinSize is true the RooDataHist
/// is filled with the functions density (function value times the
/// bin volume) rather than function value.
///
/// If showProgress is true
/// a process indicator is printed on stdout in steps of one percent,
/// which is mostly useful for the sampling of expensive functions
/// such as likelihoods

RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSet, Double_t scaleFactor,
                  Bool_t correctForBinSize, Bool_t showProgress) const
{
  // Do we have a valid histogram to use?
  if(0 == hist) {
    coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillDataHist: no valid RooDataHist to fill" << endl;
    return 0;
  }

  // Call checkObservables
  RooArgSet allDeps(*hist->get()) ;
  if (checkObservables(&allDeps)) {
    coutE(InputArguments) << "RooAbsReal::fillDataHist(" << GetName() << ") error in checkObservables, abort" << endl ;
    return hist ;
  }

  // Make deep clone of self and attach to dataset observables
  //RooArgSet* origObs = getObservables(hist) ;
  RooArgSet* cloneSet = (RooArgSet*) RooArgSet(*this).snapshot(kTRUE) ;
  RooAbsReal* theClone = (RooAbsReal*) cloneSet->find(GetName()) ;
  theClone->recursiveRedirectServers(*hist->get()) ;
  //const_cast<RooAbsReal*>(this)->recursiveRedirectServers(*hist->get()) ;

  // Iterator over all bins of RooDataHist and fill weights
  Int_t onePct = hist->numEntries()/100 ;
  if (onePct==0) {
    onePct++ ;
  }
  for (Int_t i=0 ; i<hist->numEntries() ; i++) {
    if (showProgress && (i%onePct==0)) {
      ccoutP(Eval) << "." << flush ;
    }
    const RooArgSet* obs = hist->get(i) ;
    Double_t binVal = theClone->getVal(normSet?normSet:obs)*scaleFactor ;
    if (correctForBinSize) {
      binVal*= hist->binVolume() ;
    }
    hist->set(binVal) ;
  }

  delete cloneSet ;
  //const_cast<RooAbsReal*>(this)->recursiveRedirectServers(*origObs) ;
  //delete origObs ;

  return hist;
}




////////////////////////////////////////////////////////////////////////////////
/// Create and fill a ROOT histogram TH1,TH2 or TH3 with the values of this function for the variables with given names
/// The number of bins can be controlled using the [xyz]bins parameters. For a greater degree of control
/// use the createHistogram() method below with named arguments
///
/// The caller takes ownership of the returned histogram

TH1* RooAbsReal::createHistogram(const char* varNameList, Int_t xbins, Int_t ybins, Int_t zbins) const
{
  // Parse list of variable names
  char buf[1024] ;
  strlcpy(buf,varNameList,1024) ;
  char* varName = strtok(buf,",:") ;

  RooArgSet* vars = getVariables() ;

  RooRealVar* xvar = (RooRealVar*) vars->find(varName) ;
  varName = strtok(0,",") ;
  RooRealVar* yvar = varName ? (RooRealVar*) vars->find(varName) : 0 ;
  varName = strtok(0,",") ;
  RooRealVar* zvar = varName ? (RooRealVar*) vars->find(varName) : 0 ;

  delete vars ;

  // Construct list of named arguments to pass to the implementation version of createHistogram()

  RooLinkedList argList ;
  if (xbins>0) {
    argList.Add(RooFit::Binning(xbins).Clone()) ;
  }

  if (yvar) {
    if (ybins>0) {
      argList.Add(RooFit::YVar(*yvar,RooFit::Binning(ybins)).Clone()) ;
    } else {
      argList.Add(RooFit::YVar(*yvar).Clone()) ;
    }
  }


  if (zvar) {
    if (zbins>0) {
      argList.Add(RooFit::ZVar(*zvar,RooFit::Binning(zbins)).Clone()) ;
    } else {
      argList.Add(RooFit::ZVar(*zvar).Clone()) ;
    }
  }


  // Call implementation function
  TH1* result = createHistogram(GetName(),*xvar,argList) ;

  // Delete temporary list of RooCmdArgs
  argList.Delete() ;

  return result ;
}



////////////////////////////////////////////////////////////////////////////////
/// Create and fill a ROOT histogram TH1,TH2 or TH3 with the values of this function.
///
/// This function accepts the following arguments
///
/// name -- Name of the ROOT histogram
/// xvar -- Observable to be mapped on x axis of ROOT histogram
///
/// IntrinsicBinning()                           -- Apply binning defined by function or pdf (as advertised via binBoundaries() method)
/// Binning(const char* name)                    -- Apply binning with given name to x axis of histogram
/// Binning(RooAbsBinning& binning)              -- Apply specified binning to x axis of histogram
/// Binning(int nbins, [double lo, double hi])   -- Apply specified binning to x axis of histogram
/// ConditionalObservables(const RooArgSet& set) -- Do not normalized PDF over following observables when projecting PDF into histogram
/// Scaling(Bool_t)                              -- Apply density-correction scaling (multiply by bin volume), default is kTRUE
/// Extended(Bool_t)                             -- Plot event yield instead of probability density (for extended pdfs only)
///
/// YVar(const RooAbsRealLValue& var,...)    -- Observable to be mapped on y axis of ROOT histogram
/// ZVar(const RooAbsRealLValue& var,...)    -- Observable to be mapped on z axis of ROOT histogram
///
/// The YVar() and ZVar() arguments can be supplied with optional Binning() arguments to control the binning of the Y and Z axes, e.g.
/// createHistogram("histo",x,Binning(-1,1,20), YVar(y,Binning(-1,1,30)), ZVar(z,Binning("zbinning")))
///
/// The caller takes ownership of the returned histogram

TH1 *RooAbsReal::createHistogram(const char *name, const RooAbsRealLValue& xvar,
             const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
             const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
{

  RooLinkedList l ;
  l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
  l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
  l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
  l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;

  return createHistogram(name,xvar,l) ;
}


////////////////////////////////////////////////////////////////////////////////
/// Internal method implementing createHistogram

TH1* RooAbsReal::createHistogram(const char *name, const RooAbsRealLValue& xvar, RooLinkedList& argList) const
{

  // Define configuration for this method
  RooCmdConfig pc(Form("RooAbsReal::createHistogram(%s)",GetName())) ;
  pc.defineInt("scaling","Scaling",0,1) ;
  pc.defineInt("intBinning","IntrinsicBinning",0,2) ;
  pc.defineInt("extended","Extended",0,2) ;

  pc.defineObject("compSet","SelectCompSet",0) ;
  pc.defineString("compSpec","SelectCompSpec",0) ;
  pc.defineSet("projObs","ProjectedObservables",0,0) ;
  pc.defineObject("yvar","YVar",0,0) ;
  pc.defineObject("zvar","ZVar",0,0) ;
  pc.defineMutex("SelectCompSet","SelectCompSpec") ;
  pc.defineMutex("IntrinsicBinning","Binning") ;
  pc.defineMutex("IntrinsicBinning","BinningName") ;
  pc.defineMutex("IntrinsicBinning","BinningSpec") ;
  pc.allowUndefined() ;

  // Process & check varargs
  pc.process(argList) ;
  if (!pc.ok(kTRUE)) {
    return 0 ;
  }

  RooArgList vars(xvar) ;
  RooAbsArg* yvar = static_cast<RooAbsArg*>(pc.getObject("yvar")) ;
  if (yvar) {
    vars.add(*yvar) ;
  }
  RooAbsArg* zvar = static_cast<RooAbsArg*>(pc.getObject("zvar")) ;
  if (zvar) {
    vars.add(*zvar) ;
  }

  RooArgSet* projObs = pc.getSet("projObs") ;
  RooArgSet* intObs = 0 ;

  Bool_t doScaling = pc.getInt("scaling") ;
  Int_t doIntBinning = pc.getInt("intBinning") ;
  Int_t doExtended = pc.getInt("extended") ;

  // If doExtended is two, selection is automatic, set to 1 of pdf is extended, to zero otherwise
  const RooAbsPdf* pdfSelf = dynamic_cast<const RooAbsPdf*>(this) ;
  if (!pdfSelf && doExtended>0) {
    coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-pdf object, ignored" << endl ;
    doExtended=0 ;
  }
  if (pdfSelf && doExtended==1 && pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended) {
    coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-extendable pdf, ignored" << endl ;
    doExtended=0 ;
  }
  if (pdfSelf && doExtended==2) {
    doExtended = pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended ? 0 : 1 ;
  }

  const char* compSpec = pc.getString("compSpec") ;
  const RooArgSet* compSet = (const RooArgSet*) pc.getObject("compSet") ;
  Bool_t haveCompSel = ( (compSpec && strlen(compSpec)>0) || compSet) ;

  RooBinning* intBinning(0) ;
  if (doIntBinning>0) {
    // Given RooAbsPdf* pdf and RooRealVar* obs
    list<Double_t>* bl = binBoundaries((RooRealVar&)xvar,xvar.getMin(),xvar.getMax()) ;
    if (!bl) {
      // Only emit warning when intrinsic binning is explicitly requested
      if (doIntBinning==1) {
   coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
               << ") WARNING, intrinsic model binning requested for histogram, but model does not define bin boundaries, reverting to default binning"<< endl ;
      }
    } else {
      if (doIntBinning==2) {
   coutI(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
               << ") INFO: Model has intrinsic binning definition, selecting that binning for the histogram"<< endl ;
      }
      Double_t* ba = new Double_t[bl->size()] ; int i=0 ;
      for (list<double>::iterator it=bl->begin() ; it!=bl->end() ; ++it) { ba[i++] = *it ; }
      intBinning = new RooBinning(bl->size()-1,ba) ;
      delete[] ba ;
    }
  }

  RooLinkedList argListCreate(argList) ;
  pc.stripCmdList(argListCreate,"Scaling,ProjectedObservables,IntrinsicBinning,SelectCompSet,SelectCompSpec,Extended") ;

  TH1* histo(0) ;
  if (intBinning) {
    RooCmdArg tmp = RooFit::Binning(*intBinning) ;
    argListCreate.Add(&tmp) ;
    histo = xvar.createHistogram(name,argListCreate) ;
    delete intBinning ;
  } else {
    histo = xvar.createHistogram(name,argListCreate) ;
  }

  // Do component selection here
  if (haveCompSel) {

    // Get complete set of tree branch nodes
    RooArgSet branchNodeSet ;
    branchNodeServerList(&branchNodeSet) ;

    // Discard any non-RooAbsReal nodes
    TIterator* iter = branchNodeSet.createIterator() ;
    RooAbsArg* arg ;
    while((arg=(RooAbsArg*)iter->Next())) {
      if (!dynamic_cast<RooAbsReal*>(arg)) {
   branchNodeSet.remove(*arg) ;
      }
    }
    delete iter ;

    RooArgSet* dirSelNodes ;
    if (compSet) {
      dirSelNodes = (RooArgSet*) branchNodeSet.selectCommon(*compSet) ;
    } else {
      dirSelNodes = (RooArgSet*) branchNodeSet.selectByName(compSpec) ;
    }
    if (dirSelNodes->getSize()>0) {
      coutI(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") directly selected PDF components: " << *dirSelNodes << endl ;

      // Do indirect selection and activate both
      plotOnCompSelect(dirSelNodes) ;
    } else {
      if (compSet) {
   coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection set " << *compSet << " does not match any components of p.d.f." << endl ;
      } else {
   coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection expression '" << compSpec << "' does not select any components of p.d.f." << endl ;
      }
      return 0 ;
    }
    delete dirSelNodes ;
  }

  Double_t scaleFactor(1.0) ;
  if (doExtended) {
    scaleFactor = pdfSelf->expectedEvents(vars) ;
    doScaling=kFALSE ;
  }

  fillHistogram(histo,vars,scaleFactor,intObs,doScaling,projObs,kFALSE) ;

  // Deactivate component selection
  if (haveCompSel) {
      plotOnCompSelect(0) ;
  }


  return histo ;
}


////////////////////////////////////////////////////////////////////////////////
/// Helper function for plotting of composite p.d.fs. Given
/// a set of selected components that should be plotted,
/// find all nodes that (in)directly depend on these selected
/// nodes. Mark all directly and indirecty selected nodes
/// as 'selected' using the selectComp() method

void RooAbsReal::plotOnCompSelect(RooArgSet* selNodes) const
{
  // Get complete set of tree branch nodes
  RooArgSet branchNodeSet ;
  branchNodeServerList(&branchNodeSet) ;

  // Discard any non-PDF nodes
  TIterator* iter = branchNodeSet.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (!dynamic_cast<RooAbsReal*>(arg)) {
      branchNodeSet.remove(*arg) ;
    }
  }

  // If no set is specified, restored all selection bits to kTRUE
  if (!selNodes) {
    // Reset PDF selection bits to kTRUE
    iter->Reset() ;
    while((arg=(RooAbsArg*)iter->Next())) {
      ((RooAbsReal*)arg)->selectComp(kTRUE) ;
    }
    delete iter ;
    return ;
  }


  // Add all nodes below selected nodes
  iter->Reset() ;
  TIterator* sIter = selNodes->createIterator() ;
  RooArgSet tmp ;
  while((arg=(RooAbsArg*)iter->Next())) {
    sIter->Reset() ;
    RooAbsArg* selNode ;
    while((selNode=(RooAbsArg*)sIter->Next())) {
      if (selNode->dependsOn(*arg)) {
   tmp.add(*arg,kTRUE) ;
      }
    }
  }
  delete sIter ;

  // Add all nodes that depend on selected nodes
  iter->Reset() ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (arg->dependsOn(*selNodes)) {
      tmp.add(*arg,kTRUE) ;
    }
  }

  tmp.remove(*selNodes,kTRUE) ;
  tmp.remove(*this) ;
  selNodes->add(tmp) ;
  coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") indirectly selected PDF components: " << tmp << endl ;

  // Set PDF selection bits according to selNodes
  iter->Reset() ;
  while((arg=(RooAbsArg*)iter->Next())) {
    Bool_t select = selNodes->find(arg->GetName()) ? kTRUE : kFALSE ;
    ((RooAbsReal*)arg)->selectComp(select) ;
  }

  delete iter ;
}



////////////////////////////////////////////////////////////////////////////////
/// Plot (project) PDF on specified frame. If a PDF is plotted in an empty frame, it
/// will show a unit normalized curve in the frame variable, taken at the present value
/// of other observables defined for this PDF
///
/// If a PDF is plotted in a frame in which a dataset has already been plotted, it will
/// show a projected curve integrated over all variables that were present in the shown
/// dataset except for the one on the x-axis. The normalization of the curve will also
/// be adjusted to the event count of the plotted dataset. An informational message
/// will be printed for each projection step that is performed
///
/// This function takes the following named arguments
///
/// Projection control
/// ------------------
/// Slice(const RooArgSet& set)     -- Override default projection behaviour by omittting observables listed
///                                    in set from the projection, resulting a 'slice' plot. Slicing is usually
///                                    only sensible in discrete observables. The slice is position at the 'current'
///                                    value of the observable objects
///
/// Slice(RooCategory& cat,         -- Override default projection behaviour by omittting specified category
///       const char* label)           observable from the projection, resulting in a 'slice' plot. The slice is positioned
///                                    at the given label value. Multiple Slice() commands can be given to specify slices
///                                    in multiple observables
///
/// Project(const RooArgSet& set)   -- Override default projection behaviour by projecting over observables
///                                    given in set and complete ignoring the default projection behavior. Advanced use only.
///
/// ProjWData(const RooAbsData& d)  -- Override default projection _technique_ (integration). For observables present in given dataset
///                                    projection of PDF is achieved by constructing an average over all observable values in given set.
///                                    Consult RooFit plotting tutorial for further explanation of meaning & use of this technique
///
/// ProjWData(const RooArgSet& s,   -- As above but only consider subset 's' of observables in dataset 'd' for projection through data averaging
///           const RooAbsData& d)
///
/// ProjectionRange(const char* rn) -- Override default range of projection integrals to a different range speficied by given range name.
///                                    This technique allows you to project a finite width slice in a real-valued observable
///
/// NumCPU(Int_t ncpu)              -- Number of CPUs to use simultaneously to calculate data-weighted projections (only in combination with ProjWData)
///
///
/// Misc content control
/// --------------------
/// PrintEvalErrors(Int_t numErr)   -- Control number of p.d.f evaluation errors printed per curve. A negative
///                                    value suppress output completely, a zero value will only print the error count per p.d.f component,
///                                    a positive value is will print details of each error up to numErr messages per p.d.f component.
///
/// EvalErrorValue(Double_t value)  -- Set curve points at which (pdf) evaluation error occur to specified value. By default the
///                                    function value is plotted.
///
/// Normalization(Double_t scale,   -- Adjust normalization by given scale factor. Interpretation of number depends on code: Relative:
///                ScaleType code)     relative adjustment factor, NumEvent: scale to match given number of events.
///
/// Name(const chat* name)          -- Give curve specified name in frame. Useful if curve is to be referenced later
///
/// Asymmetry(const RooCategory& c) -- Show the asymmetry of the PDF in given two-state category [F(+)-F(-)] / [F(+)+F(-)] rather than
///                                    the PDF projection. Category must have two states with indices -1 and +1 or three states with
///                                    indeces -1,0 and +1.
///
/// ShiftToZero(Bool_t flag)        -- Shift entire curve such that lowest visible point is at exactly zero. Mostly useful when
///                                    plotting -log(L) or chi^2 distributions
///
/// AddTo(const char* name,         -- Add constructed projection to already existing curve with given name and relative weight factors
///                                    double_t wgtSelf, double_t wgtOther)
///
/// Plotting control
/// ----------------
/// DrawOption(const char* opt)     -- Select ROOT draw option for resulting TGraph object
///
/// LineStyle(Int_t style)          -- Select line style by ROOT line style code, default is solid
///
/// LineColor(Int_t color)          -- Select line color by ROOT color code, default is blue
///
/// LineWidth(Int_t width)          -- Select line with in pixels, default is 3
///
/// FillStyle(Int_t style)          -- Select fill style, default is not filled. If a filled style is selected, also use VLines()
///                                    to add vertical downward lines at end of curve to ensure proper closure
/// FillColor(Int_t color)          -- Select fill color by ROOT color code
///
/// Range(const char* name)         -- Only draw curve in range defined by given name
///
/// Range(double lo, double hi)     -- Only draw curve in specified range
///
/// VLines()                        -- Add vertical lines to y=0 at end points of curve
///
/// Precision(Double_t eps)         -- Control precision of drawn curve w.r.t to scale of plot, default is 1e-3. Higher precision
///                                    will result in more and more densely spaced curve points
///
/// Invisible(Bool_t flag)           -- Add curve to frame, but do not display. Useful in combination AddTo()
///
/// VisualizeError(const RooFitResult& fitres, Double_t Z=1, Bool_t linearMethod=kTRUE)
///                                  -- Visualize the uncertainty on the parameters, as given in fitres, at 'Z' sigma'
///
/// VisualizeError(const RooFitResult& fitres, const RooArgSet& param, Double_t Z=1, Bool_t linearMethod=kTRUE) ;
///                                  -- Visualize the uncertainty on the subset of parameters 'param', as given in fitres, at 'Z' sigma'
///
///
/// Details on error band visualization
/// -----------------------------------
///
/// By default (linMethod=kTRUE) a linearized error is shown which is calculated as follows
///                                    T
/// error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
///
/// where     F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2, with f(x) the plotted curve and 'da' taken from the fit result
///       Corr(a,a') = the correlation matrix from the fit result
///                Z = requested significance 'Z sigma band'
///
/// The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters), but may
/// not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made
///
/// Alternatively (linMethod=kFALSE), a more robust error is calculated using a sampling method. In this method a number of curves
/// is calculated with variations of the parameter values, as drawn from a multi-variate Gaussian p.d.f. that is constructed
/// from the fit results covariance matrix. The error(x) is determined by calculating a central interval that capture N% of the variations
/// for each valye of x, where N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves is chosen to be such
/// that at least 30 curves are expected to be outside the N% interval, and is minimally 100 (e.g. Z=1->Ncurve=100, Z=2->Ncurve=659, Z=3->Ncurve=11111)
/// Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations, but may take (much)
/// longer to calculate

RooPlot* RooAbsReal::plotOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2,
             const RooCmdArg& arg3, const RooCmdArg& arg4,
             const RooCmdArg& arg5, const RooCmdArg& arg6,
             const RooCmdArg& arg7, const RooCmdArg& arg8,
             const RooCmdArg& arg9, const RooCmdArg& arg10) const
{
  RooLinkedList l ;
  l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
  l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
  l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
  l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;
  l.Add((TObject*)&arg9) ;  l.Add((TObject*)&arg10) ;
  return plotOn(frame,l) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Internal back-end function of plotOn() with named arguments

RooPlot* RooAbsReal::plotOn(RooPlot* frame, RooLinkedList& argList) const
{
  // Special handling here if argList contains RangeWithName argument with multiple
  // range names -- Need to translate this call into multiple calls

  RooCmdArg* rcmd = (RooCmdArg*) argList.FindObject("RangeWithName") ;
  if (rcmd && TString(rcmd->getString(0)).Contains(",")) {

    // List joint ranges as choice of normalization for all later processing
    RooCmdArg rnorm = RooFit::NormRange(rcmd->getString(0)) ;
    argList.Add(&rnorm) ;

    list<string> rlist ;

    // Separate named ranges using strtok
    char buf[1024] ;
    strlcpy(buf,rcmd->getString(0),1024) ;
    char* oneRange = strtok(buf,",") ;
    while(oneRange) {
      rlist.push_back(oneRange) ;
      oneRange = strtok(0,",") ;
    }

    for (list<string>::iterator riter=rlist.begin() ; riter!=rlist.end() ; ++riter) {
      // Process each range with a separate command with a single range to be plotted
      rcmd->setString(0,riter->c_str()) ;
      RooAbsReal::plotOn(frame,argList) ;
    }
    return frame ;

  }

  // Define configuration for this method
  RooCmdConfig pc(Form("RooAbsReal::plotOn(%s)",GetName())) ;
  pc.defineString("drawOption","DrawOption",0,"L") ;
  pc.defineString("projectionRangeName","ProjectionRange",0,"",kTRUE) ;
  pc.defineString("curveNameSuffix","CurveNameSuffix",0,"") ;
  pc.defineString("sliceCatState","SliceCat",0,"",kTRUE) ;
  pc.defineDouble("scaleFactor","Normalization",0,1.0) ;
  pc.defineObject("sliceSet","SliceVars",0) ;
  pc.defineObject("sliceCatList","SliceCat",0,0,kTRUE) ;
  pc.defineObject("projSet","Project",0) ;
  pc.defineObject("asymCat","Asymmetry",0) ;
  pc.defineDouble("precision","Precision",0,1e-3) ;
  pc.defineDouble("evalErrorVal","EvalErrorValue",0,0) ;
  pc.defineInt("doEvalError","EvalErrorValue",0,0) ;
  pc.defineInt("shiftToZero","ShiftToZero",0,0) ;
  pc.defineObject("projDataSet","ProjData",0) ;
  pc.defineObject("projData","ProjData",1) ;
  pc.defineObject("errorFR","VisualizeError",0) ;
  pc.defineDouble("errorZ","VisualizeError",0,1.) ;
  pc.defineSet("errorPars","VisualizeError",0) ;
  pc.defineInt("linearMethod","VisualizeError",0,0) ;
  pc.defineInt("binProjData","ProjData",0,0) ;
  pc.defineDouble("rangeLo","Range",0,-999.) ;
  pc.defineDouble("rangeHi","Range",1,-999.) ;
  pc.defineInt("numee","PrintEvalErrors",0,10) ;
  pc.defineInt("rangeAdjustNorm","Range",0,0) ;
  pc.defineInt("rangeWNAdjustNorm","RangeWithName",0,0) ;
  pc.defineInt("VLines","VLines",0,2) ; // 2==ExtendedWings
  pc.defineString("rangeName","RangeWithName",0,"") ;
  pc.defineString("normRangeName","NormRange",0,"") ;
  pc.defineInt("lineColor","LineColor",0,-999) ;
  pc.defineInt("lineStyle","LineStyle",0,-999) ;
  pc.defineInt("lineWidth","LineWidth",0,-999) ;
  pc.defineInt("fillColor","FillColor",0,-999) ;
  pc.defineInt("fillStyle","FillStyle",0,-999) ;
  pc.defineString("curveName","Name",0,"") ;
  pc.defineInt("curveInvisible","Invisible",0,0) ;
  pc.defineInt("showProg","ShowProgress",0,0) ;
  pc.defineInt("numCPU","NumCPU",0,1) ;
  pc.defineInt("interleave","NumCPU",1,0) ;
  pc.defineString("addToCurveName","AddTo",0,"") ;
  pc.defineDouble("addToWgtSelf","AddTo",0,1.) ;
  pc.defineDouble("addToWgtOther","AddTo",1,1.) ;
  pc.defineInt("moveToBack","MoveToBack",0,0) ;
  pc.defineMutex("SliceVars","Project") ;
  pc.defineMutex("AddTo","Asymmetry") ;
  pc.defineMutex("Range","RangeWithName") ;
  pc.defineMutex("VisualizeError","VisualizeErrorData") ;

  // Process & check varargs
  pc.process(argList) ;
  if (!pc.ok(kTRUE)) {
    return frame ;
  }

  PlotOpt o ;

  RooFitResult* errFR = (RooFitResult*) pc.getObject("errorFR") ;
  Double_t errZ = pc.getDouble("errorZ") ;
  RooArgSet* errPars = pc.getSet("errorPars") ;
  Bool_t linMethod = pc.getInt("linearMethod") ;
  if (errFR) {
    return plotOnWithErrorBand(frame,*errFR,errZ,errPars,argList,linMethod) ;
  }

  // Extract values from named arguments
  o.numee       = pc.getInt("numee") ;
  o.drawOptions = pc.getString("drawOption") ;
  o.curveNameSuffix = pc.getString("curveNameSuffix") ;
  o.scaleFactor = pc.getDouble("scaleFactor") ;
  o.projData = (const RooAbsData*) pc.getObject("projData") ;
  o.binProjData = pc.getInt("binProjData") ;
  o.projDataSet = (const RooArgSet*) pc.getObject("projDataSet") ;
  o.numCPU = pc.getInt("numCPU") ;
  o.interleave = (RooFit::MPSplit) pc.getInt("interleave") ;
  o.eeval      = pc.getDouble("evalErrorVal") ;
  o.doeeval   = pc.getInt("doEvalError") ;

  const RooArgSet* sliceSetTmp = (const RooArgSet*) pc.getObject("sliceSet") ;
  RooArgSet* sliceSet = sliceSetTmp ? ((RooArgSet*) sliceSetTmp->Clone()) : 0 ;
  const RooArgSet* projSet = (const RooArgSet*) pc.getObject("projSet") ;
  const RooAbsCategoryLValue* asymCat = (const RooAbsCategoryLValue*) pc.getObject("asymCat") ;


  // Look for category slice arguments and add them to the master slice list if found
  const char* sliceCatState = pc.getString("sliceCatState",0,kTRUE) ;
  const RooLinkedList& sliceCatList = pc.getObjectList("sliceCatList") ;
  if (sliceCatState) {

    // Make the master slice set if it doesnt exist
    if (!sliceSet) {
      sliceSet = new RooArgSet ;
    }

    // Prepare comma separated label list for parsing
    char buf[1024] ;
    strlcpy(buf,sliceCatState,1024) ;
    const char* slabel = strtok(buf,",") ;

    // Loop over all categories provided by (multiple) Slice() arguments
    TIterator* iter = sliceCatList.MakeIterator() ;
    RooCategory* scat ;
    while((scat=(RooCategory*)iter->Next())) {
      if (slabel) {
   // Set the slice position to the value indicate by slabel
   scat->setLabel(slabel) ;
   // Add the slice category to the master slice set
   sliceSet->add(*scat,kFALSE) ;
      }
      slabel = strtok(0,",") ;
    }
    delete iter ;
  }

  o.precision = pc.getDouble("precision") ;
  o.shiftToZero = (pc.getInt("shiftToZero")!=0) ;
  Int_t vlines = pc.getInt("VLines");
  if (pc.hasProcessed("Range")) {
    o.rangeLo = pc.getDouble("rangeLo") ;
    o.rangeHi = pc.getDouble("rangeHi") ;
    o.postRangeFracScale = pc.getInt("rangeAdjustNorm") ;
    if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
  } else if (pc.hasProcessed("RangeWithName")) {
    o.normRangeName = pc.getString("rangeName",0,kTRUE) ;
    o.rangeLo = frame->getPlotVar()->getMin(pc.getString("rangeName",0,kTRUE)) ;
    o.rangeHi = frame->getPlotVar()->getMax(pc.getString("rangeName",0,kTRUE)) ;
    o.postRangeFracScale = pc.getInt("rangeWNAdjustNorm") ;
    if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
  }


  // If separate normalization range was specified this overrides previous settings
  if (pc.hasProcessed("NormRange")) {
    o.normRangeName = pc.getString("normRangeName") ;
    o.postRangeFracScale = kTRUE ;
  }

  o.wmode = (vlines==2)?RooCurve::Extended:(vlines==1?RooCurve::Straight:RooCurve::NoWings) ;
  o.projectionRangeName = pc.getString("projectionRangeName",0,kTRUE) ;
  o.curveName = pc.getString("curveName",0,kTRUE) ;
  o.curveInvisible = pc.getInt("curveInvisible") ;
  o.progress = pc.getInt("showProg") ;
  o.addToCurveName = pc.getString("addToCurveName",0,kTRUE) ;
  o.addToWgtSelf = pc.getDouble("addToWgtSelf") ;
  o.addToWgtOther = pc.getDouble("addToWgtOther") ;

  if (o.addToCurveName && !frame->findObject(o.addToCurveName,RooCurve::Class())) {
    coutE(InputArguments) << "RooAbsReal::plotOn(" << GetName() << ") cannot find existing curve " << o.addToCurveName << " to add to in RooPlot" << endl ;
    return frame ;
  }

  RooArgSet projectedVars ;
  if (sliceSet) {
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have slice " << *sliceSet << endl ;

    makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,kTRUE) ;

    // Take out the sliced variables
    TIterator* iter = sliceSet->createIterator() ;
    RooAbsArg* sliceArg ;
    while((sliceArg=(RooAbsArg*)iter->Next())) {
      RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
      if (arg) {
   projectedVars.remove(*arg) ;
      } else {
   coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") slice variable "
         << sliceArg->GetName() << " was not projected anyway" << endl ;
      }
    }
    delete iter ;
  } else if (projSet) {
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have projSet " << *projSet << endl ;
    makeProjectionSet(frame->getPlotVar(),projSet,projectedVars,kFALSE) ;
  } else {
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have neither sliceSet nor projSet " << endl ;
    makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,kTRUE) ;
  }
  o.projSet = &projectedVars ;

  cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: projectedVars = " << projectedVars << endl ;


  RooPlot* ret ;
  if (!asymCat) {
    // Forward to actual calculation
    ret = RooAbsReal::plotOn(frame,o) ;
  } else {
    // Forward to actual calculation
    ret = RooAbsReal::plotAsymOn(frame,*asymCat,o) ;
  }

  delete sliceSet ;

  // Optionally adjust line/fill attributes
  Int_t lineColor = pc.getInt("lineColor") ;
  Int_t lineStyle = pc.getInt("lineStyle") ;
  Int_t lineWidth = pc.getInt("lineWidth") ;
  Int_t fillColor = pc.getInt("fillColor") ;
  Int_t fillStyle = pc.getInt("fillStyle") ;
  if (lineColor!=-999) ret->getAttLine()->SetLineColor(lineColor) ;
  if (lineStyle!=-999) ret->getAttLine()->SetLineStyle(lineStyle) ;
  if (lineWidth!=-999) ret->getAttLine()->SetLineWidth(lineWidth) ;
  if (fillColor!=-999) ret->getAttFill()->SetFillColor(fillColor) ;
  if (fillStyle!=-999) ret->getAttFill()->SetFillStyle(fillStyle) ;

  // Move last inserted object to back to drawing stack if requested
  if (pc.getInt("moveToBack") && frame->numItems()>1) {
    frame->drawBefore(frame->getObject(0)->GetName(), frame->getCurve()->GetName());
  }

  return ret ;
}




//_____________________________________________________________________________
// coverity[PASS_BY_VALUE]
RooPlot* RooAbsReal::plotOn(RooPlot *frame, PlotOpt o) const
{
  // Plotting engine function for internal use
  //
  // Plot ourselves on given frame. If frame contains a histogram, all dimensions of the plotted
  // function that occur in the previously plotted dataset are projected via partial integration,
  // otherwise no projections are performed. Optionally, certain projections can be performed
  // by summing over the values present in a provided dataset ('projData'), to correctly
  // project out data dependents that are not properly described by the PDF (e.g. per-event errors).
  //
  // The functions value can be multiplied with an optional scale factor. The interpretation
  // of the scale factor is unique for generic real functions, for PDFs there are various interpretations
  // possible, which can be selection with 'stype' (see RooAbsPdf::plotOn() for details).
  //
  // The default projection behaviour can be overriden by supplying an optional set of dependents
  // to project. For most cases, plotSliceOn() and plotProjOn() provide a more intuitive interface
  // to modify the default projection behavour.

  // Sanity checks
  if (plotSanityChecks(frame)) return frame ;

  // ProjDataVars is either all projData observables, or the user indicated subset of it
  RooArgSet projDataVars ;
  if (o.projData) {
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjData with observables = " << *o.projData->get() << endl ;
    if (o.projDataSet) {
      RooArgSet* tmp = (RooArgSet*) o.projData->get()->selectCommon(*o.projDataSet) ;
      projDataVars.add(*tmp) ;
      cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjDataSet = " << *o.projDataSet << " will only use this subset of projData" << endl ;
      delete tmp ;
    } else {
      cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") using full ProjData" << endl ;
      projDataVars.add(*o.projData->get()) ;
    }
  }

  cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") ProjDataVars = " << projDataVars << endl ;

  // Make list of variables to be projected
  RooArgSet projectedVars ;
  RooArgSet sliceSet ;
  if (o.projSet) {
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have input projSet = " << *o.projSet << endl ;
    makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,kFALSE) ;
    cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") calculated projectedVars = " << *o.projSet << endl ;

    // Print list of non-projected variables
    if (frame->getNormVars()) {
      RooArgSet *sliceSetTmp = getObservables(*frame->getNormVars()) ;

      cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") frame->getNormVars() that are also observables = " << *sliceSetTmp << endl ;

      sliceSetTmp->remove(projectedVars,kTRUE,kTRUE) ;
      sliceSetTmp->remove(*frame->getPlotVar(),kTRUE,kTRUE) ;

      if (o.projData) {
   RooArgSet* tmp = (RooArgSet*) projDataVars.selectCommon(*o.projSet) ;
   sliceSetTmp->remove(*tmp,kTRUE,kTRUE) ;
   delete tmp ;
      }

      if (sliceSetTmp->getSize()) {
   coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on "
         << frame->getPlotVar()->GetName() << " represents a slice in " << *sliceSetTmp << endl ;
      }
      sliceSet.add(*sliceSetTmp) ;
      delete sliceSetTmp ;
    }
  } else {
    makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,kTRUE) ;
  }

  cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") projectedVars = " << projectedVars << " sliceSet = " << sliceSet << endl ;


  RooArgSet* projDataNeededVars = 0 ;
  // Take out data-projected dependents from projectedVars
  if (o.projData) {
    projDataNeededVars = (RooArgSet*) projectedVars.selectCommon(projDataVars) ;
    projectedVars.remove(projDataVars,kTRUE,kTRUE) ;
  }

  // Clone the plot variable
  RooAbsReal* realVar = (RooRealVar*) frame->getPlotVar() ;
  RooArgSet* plotCloneSet = (RooArgSet*) RooArgSet(*realVar).snapshot(kTRUE) ;
  if (!plotCloneSet) {
    coutE(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Couldn't deep-clone self, abort," << endl ;
    return frame ;
  }
  RooRealVar* plotVar = (RooRealVar*) plotCloneSet->find(realVar->GetName());

  // Inform user about projections
  if (projectedVars.getSize()) {
    coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
          << " integrates over variables " << projectedVars
          << (o.projectionRangeName?Form(" in range %s",o.projectionRangeName):"") << endl;
  }
  if (projDataNeededVars && projDataNeededVars->getSize()>0) {
    coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
          << " averages using data variables " << *projDataNeededVars << endl ;
  }

  // Create projection integral
  RooArgSet* projectionCompList = 0 ;

  RooArgSet* deps = getObservables(frame->getNormVars()) ;
  deps->remove(projectedVars,kTRUE,kTRUE) ;
  if (projDataNeededVars) {
    deps->remove(*projDataNeededVars,kTRUE,kTRUE) ;
  }
  deps->remove(*plotVar,kTRUE,kTRUE) ;
  deps->add(*plotVar) ;

  // Now that we have the final set of dependents, call checkObservables()

  // WVE take out conditional observables
  if (checkObservables(deps)) {
    coutE(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") error in checkObservables, abort" << endl ;
    delete deps ;
    delete plotCloneSet ;
    if (projDataNeededVars) delete projDataNeededVars ;
    return frame ;
  }

  RooArgSet normSet(*deps) ;
  //normSet.add(projDataVars) ;

  RooAbsReal *projection = (RooAbsReal*) createPlotProjection(normSet, &projectedVars, projectionCompList, o.projectionRangeName) ;
  cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot projection object is " << projection->GetName() << endl ;
  if (dologD(Plotting)) {
    projection->printStream(ccoutD(Plotting),0,kVerbose) ;
  }

  // Always fix RooAddPdf normalizations
  RooArgSet fullNormSet(*deps) ;
  fullNormSet.add(projectedVars) ;
  if (projDataNeededVars && projDataNeededVars->getSize()>0) {
    fullNormSet.add(*projDataNeededVars) ;
  }
  RooArgSet* compSet = projection->getComponents() ;
  TIterator* iter = compSet->createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
    if (pdf) {
      pdf->selectNormalization(&fullNormSet) ;
    }
  }
  delete iter ;
  delete compSet ;


  // Apply data projection, if requested
  if (o.projData && projDataNeededVars && projDataNeededVars->getSize()>0) {

    // If data set contains more rows than needed, make reduced copy first
    RooAbsData* projDataSel = (RooAbsData*)o.projData;

    if (projDataNeededVars->getSize()<o.projData->get()->getSize()) {

      // Determine if there are any slice variables in the projection set
      RooArgSet* sliceDataSet = (RooArgSet*) sliceSet.selectCommon(*o.projData->get()) ;
      TString cutString ;
      if (sliceDataSet->getSize()>0) {
   TIterator* iter2 = sliceDataSet->createIterator() ;
   RooAbsArg* sliceVar ;
   Bool_t first(kTRUE) ;
   while((sliceVar=(RooAbsArg*)iter2->Next())) {
     if (!first) {
       cutString.Append("&&") ;
     } else {
       first=kFALSE ;
     }

     RooAbsRealLValue* real ;
     RooAbsCategoryLValue* cat ;
     if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
       cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
     } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
       cutString.Append(Form("%s==%d",cat->GetName(),cat->getIndex())) ;
     }
   }
   delete iter2 ;
      }
      delete sliceDataSet ;

      if (!cutString.IsNull()) {
   projDataSel = ((RooAbsData*)o.projData)->reduce(*projDataNeededVars,cutString) ;
   coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") reducing given projection dataset to entries with " << cutString << endl ;
      } else {
   projDataSel = ((RooAbsData*)o.projData)->reduce(*projDataNeededVars) ;
      }
      coutI(Plotting) << "RooAbsReal::plotOn(" << GetName()
            << ") only the following components of the projection data will be used: " << *projDataNeededVars << endl ;
    }

    // Request binning of unbinned projection dataset that consists exclusively of category observables
    if (!o.binProjData && dynamic_cast<RooDataSet*>(projDataSel)!=0) {

      // Determine if dataset contains only categories
      TIterator* iter2 = projDataSel->get()->createIterator() ;
      Bool_t allCat(kTRUE) ;
      RooAbsArg* arg2 ;
      while((arg2=(RooAbsArg*)iter2->Next())) {
   if (!dynamic_cast<RooCategory*>(arg2)) allCat = kFALSE ;
      }
      delete iter2 ;
      if (allCat) {
   o.binProjData = kTRUE ;
   coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") unbinned projection dataset consist only of discrete variables,"
         << " performing projection with binned copy for optimization." << endl ;

      }
    }

    // Bin projection dataset if requested
    if (o.binProjData) {
      RooAbsData* tmp = new RooDataHist(Form("%s_binned",projDataSel->GetName()),"Binned projection data",*projDataSel->get(),*projDataSel) ;
      if (projDataSel!=o.projData) delete projDataSel ;
      projDataSel = tmp ;
    }



    // Attach dataset
    projection->getVal(projDataSel->get()) ;
    projection->attachDataSet(*projDataSel) ;

    // Construct optimized data weighted average
    RooDataWeightedAverage dwa(Form("%sDataWgtAvg",GetName()),"Data Weighted average",*projection,*projDataSel,RooArgSet()/**projDataSel->get()*/,o.numCPU,o.interleave,kTRUE) ;
    //RooDataWeightedAverage dwa(Form("%sDataWgtAvg",GetName()),"Data Weighted average",*projection,*projDataSel,*projDataSel->get(),o.numCPU,o.interleave,kTRUE) ;

    // Do _not_ activate cache-and-track as necessary information to define normalization observables are not present in the underlying dataset
    dwa.constOptimizeTestStatistic(Activate,kFALSE) ;

    RooRealBinding projBind(dwa,*plotVar) ;
    RooScaledFunc scaleBind(projBind,o.scaleFactor);

    // Set default range, if not specified
    if (o.rangeLo==0 && o.rangeHi==0) {
      o.rangeLo = frame->GetXaxis()->GetXmin() ;
      o.rangeHi = frame->GetXaxis()->GetXmax() ;
    }

    // Construct name of curve for data weighed average
    TString curveName(projection->GetName()) ;
    curveName.Append(Form("_DataAvg[%s]",projDataSel->get()->contentsString().c_str())) ;
    // Append slice set specification if any
    if (sliceSet.getSize()>0) {
      curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
    }
    // Append any suffixes imported from RooAbsPdf::plotOn
    if (o.curveNameSuffix) {
      curveName.Append(o.curveNameSuffix) ;
    }

    // Curve constructor for data weighted average
    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
    RooCurve *curve = new RooCurve(projection->GetName(),projection->GetTitle(),scaleBind,
               o.rangeLo,o.rangeHi,frame->GetNbinsX(),o.precision,o.precision,o.shiftToZero,o.wmode,o.numee,o.doeeval,o.eeval) ;
    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;

    curve->SetName(curveName.Data()) ;

    // Add self to other curve if requested
    if (o.addToCurveName) {
      RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;

      // Curve constructor for sum of curves
      RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
      sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
      delete curve ;
      curve = sumCurve ;

    }

    if (o.curveName) {
      curve->SetName(o.curveName) ;
    }

    // add this new curve to the specified plot frame
    frame->addPlotable(curve, o.drawOptions, o.curveInvisible);

    if (projDataSel!=o.projData) delete projDataSel ;

  } else {

    // Set default range, if not specified
    if (o.rangeLo==0 && o.rangeHi==0) {
      o.rangeLo = frame->GetXaxis()->GetXmin() ;
      o.rangeHi = frame->GetXaxis()->GetXmax() ;
    }

    // Calculate a posteriori range fraction scaling if requested (2nd part of normalization correction for
    // result fit on subrange of data)
    if (o.postRangeFracScale) {
      if (!o.normRangeName) {
   o.normRangeName = "plotRange" ;
   plotVar->setRange("plotRange",o.rangeLo,o.rangeHi) ;
      }

      // Evaluate fractional correction integral always on full p.d.f, not component.
      Bool_t tmp = _globalSelectComp ;
      globalSelectComp(kTRUE) ;
      RooAbsReal* intFrac = projection->createIntegral(*plotVar,*plotVar,o.normRangeName) ;
      globalSelectComp(kTRUE) ;
      o.scaleFactor /= intFrac->getVal() ;
      globalSelectComp(tmp) ;
      delete intFrac ;

    }

    // create a new curve of our function using the clone to do the evaluations
    // Curve constructor for regular projections

    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
    RooCurve *curve = new RooCurve(*projection,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
               o.scaleFactor,0,o.precision,o.precision,o.shiftToZero,o.wmode,o.numee,o.doeeval,o.eeval,o.progress);
    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;



    // Set default name of curve
    TString curveName(projection->GetName()) ;
    if (sliceSet.getSize()>0) {
      curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
    }
    if (o.curveNameSuffix) {
      // Append any suffixes imported from RooAbsPdf::plotOn
      curveName.Append(o.curveNameSuffix) ;
    }
    curve->SetName(curveName.Data()) ;


    // Add self to other curve if requested
    if (o.addToCurveName) {
      RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;
      RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
      sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
      delete curve ;
      curve = sumCurve ;
    }

    // Override name of curve by user name, if specified
    if (o.curveName) {
      curve->SetName(o.curveName) ;
    }

    // add this new curve to the specified plot frame
    frame->addPlotable(curve, o.drawOptions, o.curveInvisible);
  }

  if (projDataNeededVars) delete projDataNeededVars ;
  delete deps ;
  delete projectionCompList ;
  delete plotCloneSet ;
  return frame;
}




////////////////////////////////////////////////////////////////////////////////
/// OBSOLETE -- RETAINED FOR BACKWARD COMPATIBILITY. Use the plotOn(frame,Slice(...)) instead

RooPlot* RooAbsReal::plotSliceOn(RooPlot *frame, const RooArgSet& sliceSet, Option_t* drawOptions,
             Double_t scaleFactor, ScaleType stype, const RooAbsData* projData) const
{
  RooArgSet projectedVars ;
  makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,kTRUE) ;

  // Take out the sliced variables
  TIterator* iter = sliceSet.createIterator() ;
  RooAbsArg* sliceArg ;
  while((sliceArg=(RooAbsArg*)iter->Next())) {
    RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
    if (arg) {
      projectedVars.remove(*arg) ;
    } else {
      coutI(Plotting) << "RooAbsReal::plotSliceOn(" << GetName() << ") slice variable "
            << sliceArg->GetName() << " was not projected anyway" << endl ;
    }
  }
  delete iter ;

  PlotOpt o ;
  o.drawOptions = drawOptions ;
  o.scaleFactor = scaleFactor ;
  o.stype = stype ;
  o.projData = projData ;
  o.projSet = &projectedVars ;
  return plotOn(frame,o) ;
}




//_____________________________________________________________________________
// coverity[PASS_BY_VALUE]
RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asymCat, PlotOpt o) const

{
  // Plotting engine for asymmetries. Implements the functionality if plotOn(frame,Asymmetry(...)))
  //
  // Plot asymmetry of ourselves, defined as
  //
  //   asym = f(asymCat=-1) - f(asymCat=+1) / ( f(asymCat=-1) + f(asymCat=+1) )
  //
  // on frame. If frame contains a histogram, all dimensions of the plotted
  // asymmetry function that occur in the previously plotted dataset are projected via partial integration.
  // Otherwise no projections are performed,
  //
  // The asymmetry function can be multiplied with an optional scale factor. The default projection
  // behaviour can be overriden by supplying an optional set of dependents to project.

  // Sanity checks
  if (plotSanityChecks(frame)) return frame ;

  // ProjDataVars is either all projData observables, or the user indicated subset of it
  RooArgSet projDataVars ;
  if (o.projData) {
    if (o.projDataSet) {
      RooArgSet* tmp = (RooArgSet*) o.projData->get()->selectCommon(*o.projDataSet) ;
      projDataVars.add(*tmp) ;
      delete tmp ;
    } else {
      projDataVars.add(*o.projData->get()) ;
    }
  }

  // Must depend on asymCat
  if (!dependsOn(asymCat)) {
    coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
          << ") function doesn't depend on asymmetry category " << asymCat.GetName() << endl ;
    return frame ;
  }

  // asymCat must be a signCat
  if (!asymCat.isSignType()) {
    coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
          << ") asymmetry category must have 2 or 3 states with index values -1,0,1" << endl ;
    return frame ;
  }

  // Make list of variables to be projected
  RooArgSet projectedVars ;
  RooArgSet sliceSet ;
  if (o.projSet) {
    makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,kFALSE) ;

    // Print list of non-projected variables
    if (frame->getNormVars()) {
      RooArgSet *sliceSetTmp = getObservables(*frame->getNormVars()) ;
      sliceSetTmp->remove(projectedVars,kTRUE,kTRUE) ;
      sliceSetTmp->remove(*frame->getPlotVar(),kTRUE,kTRUE) ;

      if (o.projData) {
   RooArgSet* tmp = (RooArgSet*) projDataVars.selectCommon(*o.projSet) ;
   sliceSetTmp->remove(*tmp,kTRUE,kTRUE) ;
   delete tmp ;
      }

      if (sliceSetTmp->getSize()) {
   coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on "
         << frame->getPlotVar()->GetName() << " represents a slice in " << *sliceSetTmp << endl ;
      }
      sliceSet.add(*sliceSetTmp) ;
      delete sliceSetTmp ;
    }
  } else {
    makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,kTRUE) ;
  }


  // Take out data-projected dependens from projectedVars
  RooArgSet* projDataNeededVars = 0 ;
  if (o.projData) {
    projDataNeededVars = (RooArgSet*) projectedVars.selectCommon(projDataVars) ;
    projectedVars.remove(projDataVars,kTRUE,kTRUE) ;
  }

  // Take out plotted asymmetry from projection
  if (projectedVars.find(asymCat.GetName())) {
    projectedVars.remove(*projectedVars.find(asymCat.GetName())) ;
  }

  // Clone the plot variable
  RooAbsReal* realVar = (RooRealVar*) frame->getPlotVar() ;
  RooRealVar* plotVar = (RooRealVar*) realVar->Clone() ;

  // Inform user about projections
  if (projectedVars.getSize()) {
    coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on " << plotVar->GetName()
          << " projects variables " << projectedVars << endl ;
  }
  if (projDataNeededVars && projDataNeededVars->getSize()>0) {
    coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
          << " averages using data variables "<<  *projDataNeededVars << endl ;
  }


  // Customize two copies of projection with fixed negative and positive asymmetry
  RooAbsCategoryLValue* asymPos = (RooAbsCategoryLValue*) asymCat.Clone("asym_pos") ;
  RooAbsCategoryLValue* asymNeg = (RooAbsCategoryLValue*) asymCat.Clone("asym_neg") ;
  asymPos->setIndex(1) ;
  asymNeg->setIndex(-1) ;
  RooCustomizer* custPos = new RooCustomizer(*this,"pos") ;
  RooCustomizer* custNeg = new RooCustomizer(*this,"neg") ;
  //custPos->setOwning(kTRUE) ;
  //custNeg->setOwning(kTRUE) ;
  custPos->replaceArg(asymCat,*asymPos) ;
  custNeg->replaceArg(asymCat,*asymNeg) ;
  RooAbsReal* funcPos = (RooAbsReal*) custPos->build() ;
  RooAbsReal* funcNeg = (RooAbsReal*) custNeg->build() ;

  // Create projection integral
  RooArgSet *posProjCompList, *negProjCompList ;

  // Add projDataVars to normalized dependents of projection
  // This is needed only for asymmetries (why?)
  RooArgSet depPos(*plotVar,*asymPos) ;
  RooArgSet depNeg(*plotVar,*asymNeg) ;
  depPos.add(projDataVars) ;
  depNeg.add(projDataVars) ;

  const RooAbsReal *posProj = funcPos->createPlotProjection(depPos, &projectedVars, posProjCompList, o.projectionRangeName) ;
  const RooAbsReal *negProj = funcNeg->createPlotProjection(depNeg, &projectedVars, negProjCompList, o.projectionRangeName) ;
  if (!posProj || !negProj) {
    coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") Unable to create projections, abort" << endl ;
    return frame ;
  }

  // Create a RooFormulaVar representing the asymmetry
  TString asymName(GetName()) ;
  asymName.Append("_Asym[") ;
  asymName.Append(asymCat.GetName()) ;
  asymName.Append("]") ;
  TString asymTitle(asymCat.GetName()) ;
  asymTitle.Append(" Asymmetry of ") ;
  asymTitle.Append(GetTitle()) ;
  RooFormulaVar* funcAsym = new RooFormulaVar(asymName,asymTitle,"(@0-@1)/(@0+@1)",RooArgSet(*posProj,*negProj)) ;

  if (o.projData) {

    // If data set contains more rows than needed, make reduced copy first
    RooAbsData* projDataSel = (RooAbsData*)o.projData;
    if (projDataNeededVars && projDataNeededVars->getSize()<o.projData->get()->getSize()) {

      // Determine if there are any slice variables in the projection set
      RooArgSet* sliceDataSet = (RooArgSet*) sliceSet.selectCommon(*o.projData->get()) ;
      TString cutString ;
      if (sliceDataSet->getSize()>0) {
   TIterator* iter = sliceDataSet->createIterator() ;
   RooAbsArg* sliceVar ;
   Bool_t first(kTRUE) ;
   while((sliceVar=(RooAbsArg*)iter->Next())) {
     if (!first) {
       cutString.Append("&&") ;
     } else {
       first=kFALSE ;
     }

     RooAbsRealLValue* real ;
     RooAbsCategoryLValue* cat ;
     if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
       cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
     } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
       cutString.Append(Form("%s==%d",cat->GetName(),cat->getIndex())) ;
     }
   }
   delete iter ;
      }
      delete sliceDataSet ;

      if (!cutString.IsNull()) {
   projDataSel = ((RooAbsData*)o.projData)->reduce(*projDataNeededVars,cutString) ;
   coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
         << ") reducing given projection dataset to entries with " << cutString << endl ;
      } else {
   projDataSel = ((RooAbsData*)o.projData)->reduce(*projDataNeededVars) ;
      }
      coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
            << ") only the following components of the projection data will be used: " << *projDataNeededVars << endl ;
    }


    RooDataWeightedAverage dwa(Form("%sDataWgtAvg",GetName()),"Data Weighted average",*funcAsym,*projDataSel,RooArgSet()/**projDataSel->get()*/,o.numCPU,o.interleave,kTRUE) ;
    //RooDataWeightedAverage dwa(Form("%sDataWgtAvg",GetName()),"Data Weighted average",*funcAsym,*projDataSel,*projDataSel->get(),o.numCPU,o.interleave,kTRUE) ;
    dwa.constOptimizeTestStatistic(Activate) ;

    RooRealBinding projBind(dwa,*plotVar) ;

    ((RooAbsReal*)posProj)->attachDataSet(*projDataSel) ;
    ((RooAbsReal*)negProj)->attachDataSet(*projDataSel) ;

    RooScaledFunc scaleBind(projBind,o.scaleFactor);

    // Set default range, if not specified
    if (o.rangeLo==0 && o.rangeHi==0) {
      o.rangeLo = frame->GetXaxis()->GetXmin() ;
      o.rangeHi = frame->GetXaxis()->GetXmax() ;
    }

    // Construct name of curve for data weighed average
    TString curveName(funcAsym->GetName()) ;
    curveName.Append(Form("_DataAvg[%s]",projDataSel->get()->contentsString().c_str())) ;
    // Append slice set specification if any
    if (sliceSet.getSize()>0) {
      curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
    }
    // Append any suffixes imported from RooAbsPdf::plotOn
    if (o.curveNameSuffix) {
      curveName.Append(o.curveNameSuffix) ;
    }


    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
    RooCurve *curve = new RooCurve(funcAsym->GetName(),funcAsym->GetTitle(),scaleBind,
               o.rangeLo,o.rangeHi,frame->GetNbinsX(),o.precision,o.precision,kFALSE,o.wmode,o.numee,o.doeeval,o.eeval) ;
    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;

    dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;
    // add this new curve to the specified plot frame
    frame->addPlotable(curve, o.drawOptions);

    ccoutW(Eval) << endl ;

    if (projDataSel!=o.projData) delete projDataSel ;

  } else {

    // Set default range, if not specified
    if (o.rangeLo==0 && o.rangeHi==0) {
      o.rangeLo = frame->GetXaxis()->GetXmin() ;
      o.rangeHi = frame->GetXaxis()->GetXmax() ;
    }

    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
    RooCurve* curve= new RooCurve(*funcAsym,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
              o.scaleFactor,0,o.precision,o.precision,kFALSE,o.wmode,o.numee,o.doeeval,o.eeval);
    RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;

    dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;


    // Set default name of curve
    TString curveName(funcAsym->GetName()) ;
    if (sliceSet.getSize()>0) {
      curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
    }
    if (o.curveNameSuffix) {
      // Append any suffixes imported from RooAbsPdf::plotOn
      curveName.Append(o.curveNameSuffix) ;
    }
    curve->SetName(curveName.Data()) ;

    // add this new curve to the specified plot frame
    frame->addPlotable(curve, o.drawOptions);

  }

  // Cleanup
  delete custPos ;
  delete custNeg ;
  delete funcPos ;
  delete funcNeg ;
  delete posProjCompList ;
  delete negProjCompList ;
  delete asymPos ;
  delete asymNeg ;
  delete funcAsym ;

  delete plotVar ;

  return frame;
}



////////////////////////////////////////////////////////////////////////////////
/// Calculate error on self by propagated errors on parameters with correlations as given by fit result
/// The linearly propagated error is calculated as follows
///                                    T
/// error(x) = F_a(x) * Corr(a,a') F_a'(x)
///
/// where     F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2, with f(x) this function and 'da' taken from the fit result
///       Corr(a,a') = the correlation matrix from the fit result
///

Double_t RooAbsReal::getPropagatedError(const RooFitResult& fr)
{

  // Clone self for internal use
  RooAbsReal* cloneFunc = (RooAbsReal*) cloneTree() ;
  RooArgSet* errorParams = cloneFunc->getObservables(fr.floatParsFinal()) ;
  RooArgSet* nset = cloneFunc->getParameters(*errorParams) ;

  // Make list of parameter instances of cloneFunc in order of error matrix
  RooArgList paramList ;
  const RooArgList& fpf = fr.floatParsFinal() ;
  vector<int> fpf_idx ;
  for (Int_t i=0 ; i<fpf.getSize() ; i++) {
    RooAbsArg* par = errorParams->find(fpf[i].GetName()) ;
    if (par) {
      paramList.add(*par) ;
      fpf_idx.push_back(i) ;
    }
  }

  vector<Double_t> plusVar, minusVar ;

  // Create vector of plus,minus variations for each parameter
  TMatrixDSym V(paramList.getSize()==fr.floatParsFinal().getSize()?
      fr.covarianceMatrix():
      fr.reducedCovarianceMatrix(paramList)) ;

  for (Int_t ivar=0 ; ivar<paramList.getSize() ; ivar++) {

    RooRealVar& rrv = (RooRealVar&)fpf[fpf_idx[ivar]] ;

    Double_t cenVal = rrv.getVal() ;
    Double_t errVal = sqrt(V(ivar,ivar)) ;

    // Make Plus variation
    ((RooRealVar*)paramList.at(ivar))->setVal(cenVal+errVal) ;
    plusVar.push_back(cloneFunc->getVal(nset)) ;

    // Make Minus variation
    ((RooRealVar*)paramList.at(ivar))->setVal(cenVal-errVal) ;
    minusVar.push_back(cloneFunc->getVal(nset)) ;

    ((RooRealVar*)paramList.at(ivar))->setVal(cenVal) ;
  }

  TMatrixDSym C(paramList.getSize()) ;
  vector<double> errVec(paramList.getSize()) ;
  for (int i=0 ; i<paramList.getSize() ; i++) {
    errVec[i] = sqrt(V(i,i)) ;
    for (int j=i ; j<paramList.getSize() ; j++) {
      C(i,j) = V(i,j)/sqrt(V(i,i)*V(j,j)) ;
      C(j,i) = C(i,j) ;
    }
  }

  // Make vector of variations
  TVectorD F(plusVar.size()) ;
  for (unsigned int j=0 ; j<plusVar.size() ; j++) {
    F[j] = (plusVar[j]-minusVar[j])/2 ;
  }

  // Calculate error in linear approximation from variations and correlation coefficient
  Double_t sum = F*(C*F) ;

  delete cloneFunc ;
  delete errorParams ;
  delete nset ;

  return sqrt(sum) ;
}










////////////////////////////////////////////////////////////////////////////////
/// Plot function or p.d.f. on frame with support for visualization of the uncertainty encoded in the given fit result fr.
/// If params is non-zero, only the subset of the parameters in fr that occur in params is considered for the error evaluation
/// Argument argList can contain any RooCmdArg named argument that can be applied to a regular plotOn() operation
///
/// By default (linMethod=kTRUE) a linearized error is shown which is calculated as follows
///                                    T
/// error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
///
/// where     F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2, with f(x) the plotted curve and 'da' taken from the fit result
///       Corr(a,a') = the correlation matrix from the fit result
///                Z = requested signifance 'Z sigma band'
///
/// The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters), but may
/// not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made
///
/// Alternatively, a more robust error is calculated using a sampling method. In this method a number of curves
/// is calculated with variations of the parameter values, as drawn from a multi-variate Gaussian p.d.f. that is constructed
/// from the fit results covariance matrix. The error(x) is determined by calculating a central interval that capture N% of the variations
/// for each valye of x, where N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves is chosen to be such
/// that at least 30 curves are expected to be outside the N% interval, and is minimally 100 (e.g. Z=1->Ncurve=100, Z=2->Ncurve=659, Z=3->Ncurve=11111)
/// Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations

RooPlot* RooAbsReal::plotOnWithErrorBand(RooPlot* frame,const RooFitResult& fr, Double_t Z,const RooArgSet* params, const RooLinkedList& argList, Bool_t linMethod) const
{
  RooLinkedList plotArgListTmp(argList) ;
  RooCmdConfig pc(Form("RooAbsPdf::plotOn(%s)",GetName())) ;
  pc.stripCmdList(plotArgListTmp,"VisualizeError,MoveToBack") ;

  // Strip any 'internal normalization' arguments from list
  RooLinkedList plotArgList ;
  RooFIter iter = plotArgListTmp.fwdIterator() ;
  RooCmdArg* cmd ;
  while ((cmd=(RooCmdArg*)iter.next())) {
    if (std::string("Normalization")==cmd->GetName()) {
      if (((RooCmdArg*)cmd)->getInt(1)!=0) {
      } else {
   plotArgList.Add(cmd) ;
      }
    } else {
      plotArgList.Add(cmd) ;
    }
  }

  // Generate central value curve
  RooLinkedList tmp(plotArgList) ;
  plotOn(frame,tmp) ;
  RooCurve* cenCurve = frame->getCurve() ;
  frame->remove(0,kFALSE) ;

  RooCurve* band(0) ;
  if (!linMethod) {

    // *** Interval method ***
    //
    // Make N variations of parameters samples from V and visualize N% central interval where N% is defined from Z

    // Clone self for internal use
    RooAbsReal* cloneFunc = (RooAbsReal*) cloneTree() ;
    RooArgSet* cloneParams = cloneFunc->getObservables(fr.floatParsFinal()) ;
    RooArgSet* errorParams = params?((RooArgSet*)cloneParams->selectCommon(*params)):cloneParams ;

    // Generate 100 random parameter points distributed according to fit result covariance matrix
    RooAbsPdf* paramPdf = fr.createHessePdf(*errorParams) ;
    Int_t n = Int_t(100./TMath::Erfc(Z/sqrt(2.))) ;
    if (n<100) n=100 ;

    coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") INFO: visualizing " << Z << "-sigma uncertainties in parameters "
        << *errorParams << " from fit result " << fr.GetName() << " using " << n << " samplings." << endl ;

    // Generate variation curves with above set of parameter values
    Double_t ymin = frame->GetMinimum() ;
    Double_t ymax = frame->GetMaximum() ;
    RooDataSet* d = paramPdf->generate(*errorParams,n) ;
    vector<RooCurve*> cvec ;
    for (int i=0 ; i<d->numEntries() ; i++) {
      *cloneParams = (*d->get(i)) ;
      RooLinkedList tmp2(plotArgList) ;
      cloneFunc->plotOn(frame,tmp2) ;
      cvec.push_back(frame->getCurve()) ;
      frame->remove(0,kFALSE) ;
    }
    frame->SetMinimum(ymin) ;
    frame->SetMaximum(ymax) ;


    // Generate upper and lower curve points from 68% interval around each point of central curve
    band = cenCurve->makeErrorBand(cvec,Z) ;

    // Cleanup
    delete paramPdf ;
    delete cloneFunc ;
    for (vector<RooCurve*>::iterator i=cvec.begin() ; i!=cvec.end() ; i++) {
      delete (*i) ;
    }

  } else {

    // *** Linear Method ***
    //
    // Make a one-sigma up- and down fluctation for each parameter and visualize
    // a from a linearized calculation as follows
    //
    //   error(x) = F(a) C_aa' F(a')
    //
    //   Where F(a) = (f(x,a+da) - f(x,a-da))/2
    //   and C_aa' is the correlation matrix

    // Clone self for internal use
    RooAbsReal* cloneFunc = (RooAbsReal*) cloneTree() ;
    RooArgSet* cloneParams = cloneFunc->getObservables(fr.floatParsFinal()) ;
    RooArgSet* errorParams = params?((RooArgSet*)cloneParams->selectCommon(*params)):cloneParams ;


    // Make list of parameter instances of cloneFunc in order of error matrix
    RooArgList paramList ;
    const RooArgList& fpf = fr.floatParsFinal() ;
    vector<int> fpf_idx ;
    for (Int_t i=0 ; i<fpf.getSize() ; i++) {
      RooAbsArg* par = errorParams->find(fpf[i].GetName()) ;
      if (par) {
   paramList.add(*par) ;
   fpf_idx.push_back(i) ;
      }
    }

    vector<RooCurve*> plusVar, minusVar ;

    // Create vector of plus,minus variations for each parameter

    TMatrixDSym V(paramList.getSize()==fr.floatParsFinal().getSize()?
        fr.covarianceMatrix():
        fr.reducedCovarianceMatrix(paramList)) ;


    for (Int_t ivar=0 ; ivar<paramList.getSize() ; ivar++) {

      RooRealVar& rrv = (RooRealVar&)fpf[fpf_idx[ivar]] ;

      Double_t cenVal = rrv.getVal() ;
      Double_t errVal = sqrt(V(ivar,ivar)) ;

      // Make Plus variation
      ((RooRealVar*)paramList.at(ivar))->setVal(cenVal+Z*errVal) ;


      RooLinkedList tmp2(plotArgList) ;
      cloneFunc->plotOn(frame,tmp2) ;
      plusVar.push_back(frame->getCurve()) ;
      frame->remove(0,kFALSE) ;


      // Make Minus variation
      ((RooRealVar*)paramList.at(ivar))->setVal(cenVal-Z*errVal) ;
      RooLinkedList tmp3(plotArgList) ;
      cloneFunc->plotOn(frame,tmp3) ;
      minusVar.push_back(frame->getCurve()) ;
      frame->remove(0,kFALSE) ;

      ((RooRealVar*)paramList.at(ivar))->setVal(cenVal) ;
    }

    TMatrixDSym C(paramList.getSize()) ;
    vector<double> errVec(paramList.getSize()) ;
    for (int i=0 ; i<paramList.getSize() ; i++) {
      errVec[i] = sqrt(V(i,i)) ;
      for (int j=i ; j<paramList.getSize() ; j++) {
   C(i,j) = V(i,j)/sqrt(V(i,i)*V(j,j)) ;
   C(j,i) = C(i,j) ;
      }
    }

    band = cenCurve->makeErrorBand(plusVar,minusVar,C,Z) ;


    // Cleanup
    delete cloneFunc ;
    for (vector<RooCurve*>::iterator i=plusVar.begin() ; i!=plusVar.end() ; i++) {
      delete (*i) ;
    }
    for (vector<RooCurve*>::iterator i=minusVar.begin() ; i!=minusVar.end() ; i++) {
      delete (*i) ;
    }

  }

  delete cenCurve ;
  if (!band) return frame ;

  // Define configuration for this method
  pc.defineString("drawOption","DrawOption",0,"F") ;
  pc.defineString("curveNameSuffix","CurveNameSuffix",0,"") ;
  pc.defineInt("lineColor","LineColor",0,-999) ;
  pc.defineInt("lineStyle","LineStyle",0,-999) ;
  pc.defineInt("lineWidth","LineWidth",0,-999) ;
  pc.defineInt("fillColor","FillColor",0,-999) ;
  pc.defineInt("fillStyle","FillStyle",0,-999) ;
  pc.defineString("curveName","Name",0,"") ;
  pc.defineInt("curveInvisible","Invisible",0,0) ;
  pc.defineInt("moveToBack","MoveToBack",0,0) ;
  pc.allowUndefined() ;

  // Process & check varargs
  pc.process(argList) ;
  if (!pc.ok(kTRUE)) {
    return frame ;
  }

  // Insert error band in plot frame
  frame->addPlotable(band,pc.getString("drawOption"),pc.getInt("curveInvisible")) ;


  // Optionally adjust line/fill attributes
  Int_t lineColor = pc.getInt("lineColor") ;
  Int_t lineStyle = pc.getInt("lineStyle") ;
  Int_t lineWidth = pc.getInt("lineWidth") ;
  Int_t fillColor = pc.getInt("fillColor") ;
  Int_t fillStyle = pc.getInt("fillStyle") ;
  if (lineColor!=-999) frame->getAttLine()->SetLineColor(lineColor) ;
  if (lineStyle!=-999) frame->getAttLine()->SetLineStyle(lineStyle) ;
  if (lineWidth!=-999) frame->getAttLine()->SetLineWidth(lineWidth) ;
  if (fillColor!=-999) frame->getAttFill()->SetFillColor(fillColor) ;
  if (fillStyle!=-999) frame->getAttFill()->SetFillStyle(fillStyle) ;

  // Adjust name if requested
  if (pc.getString("curveName",0,kTRUE)) {
    band->SetName(pc.getString("curveName",0,kTRUE)) ;
  } else if (pc.getString("curveNameSuffix",0,kTRUE)) {
    TString name(band->GetName()) ;
    name.Append(pc.getString("curveNameSuffix",0,kTRUE)) ;
    band->SetName(name.Data()) ;
  }

  // Move last inserted object to back to drawing stack if requested
  if (pc.getInt("moveToBack") && frame->numItems()>1) {
    frame->drawBefore(frame->getObject(0)->GetName(), frame->getCurve()->GetName());
  }


  return frame ;
}




////////////////////////////////////////////////////////////////////////////////
/// Utility function for plotOn(), perform general sanity check on frame to ensure safe plotting operations

Bool_t RooAbsReal::plotSanityChecks(RooPlot* frame) const
{
  // check that we are passed a valid plot frame to use
  if(0 == frame) {
    coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
    return kTRUE;
  }

  // check that this frame knows what variable to plot
  RooAbsReal* var = frame->getPlotVar() ;
  if(!var) {
    coutE(Plotting) << ClassName() << "::" << GetName()
    << ":plotOn: frame does not specify a plot variable" << endl;
    return kTRUE;
  }

  // check that the plot variable is not derived
  if(!dynamic_cast<RooAbsRealLValue*>(var)) {
    coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: cannot plot variable \""
          << var->GetName() << "\" of type " << var->ClassName() << endl;
    return kTRUE;
  }

  // check if we actually depend on the plot variable
  if(!this->dependsOn(*var)) {
    coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: WARNING: variable is not an explicit dependent: "
          << var->GetName() << endl;
  }

  return kFALSE ;
}




////////////////////////////////////////////////////////////////////////////////
/// Utility function for plotOn() that constructs the set of
/// observables to project when plotting ourselves as function of
/// 'plotVar'. 'allVars' is the list of variables that must be
/// projected, but may contain variables that we do not depend on. If
/// 'silent' is cleared, warnings about inconsistent input parameters
/// will be printed.

void RooAbsReal::makeProjectionSet(const RooAbsArg* plotVar, const RooArgSet* allVars,
               RooArgSet& projectedVars, Bool_t silent) const
{
  cxcoutD(Plotting) << "RooAbsReal::makeProjectionSet(" << GetName() << ") plotVar = " << plotVar->GetName()
          << " allVars = " << (allVars?(*allVars):RooArgSet()) << endl ;

  projectedVars.removeAll() ;
  if (!allVars) return ;

  // Start out with suggested list of variables
  projectedVars.add(*allVars) ;

  // Take out plot variable
  RooAbsArg *found= projectedVars.find(plotVar->GetName());
  if(found) {
    projectedVars.remove(*found);

    // Take out eventual servers of plotVar
    RooArgSet* plotServers = plotVar->getObservables(&projectedVars) ;
    TIterator* psIter = plotServers->createIterator() ;
    RooAbsArg* ps ;
    while((ps=(RooAbsArg*)psIter->Next())) {
      RooAbsArg* tmp = projectedVars.find(ps->GetName()) ;
      if (tmp) {
   cxcoutD(Plotting) << "RooAbsReal::makeProjectionSet(" << GetName() << ") removing " << tmp->GetName()
           << " from projection set because it a server of " << plotVar->GetName() << endl ;
   projectedVars.remove(*tmp) ;
      }
    }
    delete psIter ;
    delete plotServers ;

    if (!silent) {
      coutW(Plotting) << "RooAbsReal::plotOn(" << GetName()
            << ") WARNING: cannot project out frame variable ("
            << found->GetName() << "), ignoring" << endl ;
    }
  }

  // Take out all non-dependents of function
  TIterator* iter = allVars->createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (!dependsOnValue(*arg)) {
      projectedVars.remove(*arg,kTRUE) ;

      cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName()
         << ") function doesn't depend on projection variable "
         << arg->GetName() << ", ignoring" << endl ;
    }
  }
  delete iter ;
}




////////////////////////////////////////////////////////////////////////////////
/// If true, the current pdf is a selected component (for use in plotting)

Bool_t RooAbsReal::isSelectedComp() const
{
  return _selectComp || _globalSelectComp ;
}



////////////////////////////////////////////////////////////////////////////////
/// Global switch controlling the activation of the selectComp() functionality

void RooAbsReal::globalSelectComp(Bool_t flag)
{
  _globalSelectComp = flag ;
}




////////////////////////////////////////////////////////////////////////////////
/// Create an interface adaptor f(vars) that binds us to the specified variables
/// (in arbitrary order). For example, calling bindVars({x1,x3}) on an object
/// F(x1,x2,x3,x4) returns an object f(x1,x3) that is evaluated using the
/// current values of x2 and x4. The caller takes ownership of the returned adaptor.

RooAbsFunc *RooAbsReal::bindVars(const RooArgSet &vars, const RooArgSet* nset, Bool_t clipInvalid) const
{
  RooAbsFunc *binding= new RooRealBinding(*this,vars,nset,clipInvalid);
  if(binding && !binding->isValid()) {
    coutE(InputArguments) << ClassName() << "::" << GetName() << ":bindVars: cannot bind to " << vars << endl ;
    delete binding;
    binding= 0;
  }
  return binding;
}



////////////////////////////////////////////////////////////////////////////////
/// Copy the cached value of another RooAbsArg to our cache.
/// Warning: This function copies the cached values of source,
/// it is the callers responsibility to make sure the cache is clean

void RooAbsReal::copyCache(const RooAbsArg* source, Bool_t /*valueOnly*/, Bool_t setValDirty)
{
  RooAbsReal* other = static_cast<RooAbsReal*>(const_cast<RooAbsArg*>(source)) ;

  if (!other->_treeVar) {
    _value = other->_value ;
  } else {
    if (source->getAttribute("FLOAT_TREE_BRANCH")) {
      _value = other->_floatValue ;
    } else if (source->getAttribute("INTEGER_TREE_BRANCH")) {
      _value = other->_intValue ;
    } else if (source->getAttribute("BYTE_TREE_BRANCH")) {
      _value = other->_byteValue ;
    } else if (source->getAttribute("BOOL_TREE_BRANCH")) {
      _value = other->_boolValue ;
    } else if (source->getAttribute("SIGNEDBYTE_TREE_BRANCH")) {
      _value = other->_sbyteValue ;
    } else if (source->getAttribute("UNSIGNED_INTEGER_TREE_BRANCH")) {
      _value = other->_uintValue ;
    }
  }
  if (setValDirty) {
    setValueDirty() ;
  }
}



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

void RooAbsReal::attachToVStore(RooVectorDataStore& vstore)
{
  RooVectorDataStore::RealVector* rv = vstore.addReal(this) ;
  rv->setBuffer(this,&_value) ;
}




////////////////////////////////////////////////////////////////////////////////
/// Attach object to a branch of given TTree. By default it will
/// register the internal value cache RooAbsReal::_value as branch
/// buffer for a Double_t tree branch with the same name as this
/// object. If no Double_t branch is found with the name of this
/// object, this method looks for a Float_t Int_t, UChar_t and UInt_t
/// branch in that order. If any of these are found the buffer for
/// that branch is set to a correctly typed conversion buffer in this
/// RooRealVar.  A flag is set that will cause copyCache to copy the
/// object value from the appropriate conversion buffer instead of
/// the _value buffer.

void RooAbsReal::attachToTree(TTree& t, Int_t bufSize)
{
  // First determine if branch is taken
  TString cleanName(cleanBranchName()) ;
  TBranch* branch = t.GetBranch(cleanName) ;
  if (branch) {

    // Determine if existing branch is Float_t or Double_t
    TLeaf* leaf = (TLeaf*)branch->GetListOfLeaves()->At(0) ;

    // Check that leaf is _not_ an array
    Int_t dummy ;
    TLeaf* counterLeaf = leaf->GetLeafCounter(dummy) ;
    if (counterLeaf) {
      coutE(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") ERROR: TTree branch " << GetName()
        << " is an array and cannot be attached to a RooAbsReal" << endl ;
      return ;
    }

    TString typeName(leaf->GetTypeName()) ;

    if (!typeName.CompareTo("Float_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree Float_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("FLOAT_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_floatValue) ;
    } else if (!typeName.CompareTo("Int_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree Int_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("INTEGER_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_intValue) ;
    } else if (!typeName.CompareTo("UChar_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree UChar_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("BYTE_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_byteValue) ;
    } else if (!typeName.CompareTo("Bool_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree Bool_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("BOOL_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_boolValue) ;
    } else if (!typeName.CompareTo("Char_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree Char_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("SIGNEDBYTE_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_sbyteValue) ;
    }  else if (!typeName.CompareTo("UInt_t")) {
      coutI(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") TTree UInt_t branch " << GetName()
        << " will be converted to double precision" << endl ;
      setAttribute("UNSIGNED_INTEGER_TREE_BRANCH",kTRUE) ;
      _treeVar = kTRUE ;
      t.SetBranchAddress(cleanName,&_uintValue) ;
    } else if (!typeName.CompareTo("Double_t")) {
      t.SetBranchAddress(cleanName,&_value) ;
    } else {
      coutE(InputArguments) << "RooAbsReal::attachToTree(" << GetName() << ") data type " << typeName << " is not supported" << endl ;
    }

    if (branch->GetCompressionLevel()<0) {
      // cout << "RooAbsReal::attachToTree(" << GetName() << ") Fixing compression level of branch " << cleanName << endl ;
      branch->SetCompressionLevel(1) ;
    }

//      cout << "RooAbsReal::attachToTree(" << cleanName << "): branch already exists in tree " << (void*)&t << ", changing address" << endl ;

  } else {

    TString format(cleanName);
    format.Append("/D");
    branch = t.Branch(cleanName, &_value, (const Text_t*)format, bufSize);
    branch->SetCompressionLevel(1) ;
    //      cout << "RooAbsReal::attachToTree(" << cleanName << "): creating new branch in tree " << (void*)&t << endl ;
  }

}



////////////////////////////////////////////////////////////////////////////////
/// Fill the tree branch that associated with this object with its current value

void RooAbsReal::fillTreeBranch(TTree& t)
{
  // First determine if branch is taken
  TBranch* branch = t.GetBranch(cleanBranchName()) ;
  if (!branch) {
    coutE(Eval) << "RooAbsReal::fillTreeBranch(" << GetName() << ") ERROR: not attached to tree: " << cleanBranchName() << endl ;
    assert(0) ;
  }
  branch->Fill() ;

}



////////////////////////////////////////////////////////////////////////////////
/// (De)Activate associated tree branch

void RooAbsReal::setTreeBranchStatus(TTree& t, Bool_t active)
{
  TBranch* branch = t.GetBranch(cleanBranchName()) ;
  if (branch) {
    t.SetBranchStatus(cleanBranchName(),active?1:0) ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Create a RooRealVar fundamental object with our properties. The new
/// object will be created without any fit limits.

RooAbsArg *RooAbsReal::createFundamental(const char* newname) const
{
  RooRealVar *fund= new RooRealVar(newname?newname:GetName(),GetTitle(),_value,getUnit());
  fund->removeRange();
  fund->setPlotLabel(getPlotLabel());
  fund->setAttribute("fundamentalCopy");
  return fund;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for use in getAnalyticalIntegral(). If the
/// content of proxy 'a' occurs in set 'allDeps' then the argument
/// held in 'a' is copied from allDeps to analDeps

Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
               const RooArgProxy& a) const
{
  TList nameList ;
  nameList.Add(new TObjString(a.absArg()->GetName())) ;
  Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
  nameList.Delete() ;
  return result ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for use in getAnalyticalIntegral(). If the
/// contents of proxies a,b occur in set 'allDeps' then the arguments
/// held in a,b are copied from allDeps to analDeps

Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
               const RooArgProxy& a, const RooArgProxy& b) const
{
  TList nameList ;
  nameList.Add(new TObjString(a.absArg()->GetName())) ;
  nameList.Add(new TObjString(b.absArg()->GetName())) ;
  Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
  nameList.Delete() ;
  return result ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for use in getAnalyticalIntegral(). If the
/// contents of proxies a,b,c occur in set 'allDeps' then the arguments
/// held in a,b,c are copied from allDeps to analDeps

Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
               const RooArgProxy& a, const RooArgProxy& b,
               const RooArgProxy& c) const
{
  TList nameList ;
  nameList.Add(new TObjString(a.absArg()->GetName())) ;
  nameList.Add(new TObjString(b.absArg()->GetName())) ;
  nameList.Add(new TObjString(c.absArg()->GetName())) ;
  Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
  nameList.Delete() ;
  return result ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for use in getAnalyticalIntegral(). If the
/// contents of proxies a,b,c,d occur in set 'allDeps' then the arguments
/// held in a,b,c,d are copied from allDeps to analDeps

Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
               const RooArgProxy& a, const RooArgProxy& b,
               const RooArgProxy& c, const RooArgProxy& d) const
{
  TList nameList ;
  nameList.Add(new TObjString(a.absArg()->GetName())) ;
  nameList.Add(new TObjString(b.absArg()->GetName())) ;
  nameList.Add(new TObjString(c.absArg()->GetName())) ;
  nameList.Add(new TObjString(d.absArg()->GetName())) ;
  Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
  nameList.Delete() ;
  return result ;
}


////////////////////////////////////////////////////////////////////////////////
/// Utility function for use in getAnalyticalIntegral(). If the
/// contents of 'refset' occur in set 'allDeps' then the arguments
/// held in 'refset' are copied from allDeps to analDeps.

Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
              const RooArgSet& refset) const
{
  TList nameList ;
  TIterator* iter = refset.createIterator() ;
  RooAbsArg* arg ;
  while ((arg=(RooAbsArg*)iter->Next())) {
    nameList.Add(new TObjString(arg->GetName())) ;
  }
  delete iter ;

  Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
  nameList.Delete() ;
  return result ;
}



////////////////////////////////////////////////////////////////////////////////
/// Check if allArgs contains matching elements for each name in nameList. If it does,
/// add the corresponding args from allArgs to matchedArgs and return kTRUE. Otherwise
/// return kFALSE and do not change matchedArgs.

Bool_t RooAbsReal::matchArgsByName(const RooArgSet &allArgs, RooArgSet &matchedArgs,
              const TList &nameList) const
{
  RooArgSet matched("matched");
  TIterator *iterator= nameList.MakeIterator();
  TObjString *name = 0;
  Bool_t isMatched(kTRUE);
  while((isMatched && (name= (TObjString*)iterator->Next()))) {
    RooAbsArg *found= allArgs.find(name->String().Data());
    if(found) {
      matched.add(*found);
    }
    else {
      isMatched= kFALSE;
    }
  }

  // nameList may not contain multiple entries with the same name
  // that are both matched
  if (isMatched && (matched.getSize()!=nameList.GetSize())) {
    isMatched = kFALSE ;
  }

  delete iterator;
  if(isMatched) matchedArgs.add(matched);
  return isMatched;
}



////////////////////////////////////////////////////////////////////////////////
/// Returns the default numeric integration configuration for all RooAbsReals

RooNumIntConfig* RooAbsReal::defaultIntegratorConfig()
{
  return &RooNumIntConfig::defaultConfig() ;
}


////////////////////////////////////////////////////////////////////////////////
/// Returns the specialized integrator configuration for _this_ RooAbsReal.
/// If this object has no specialized configuration, a null pointer is returned.

RooNumIntConfig* RooAbsReal::specialIntegratorConfig() const
{
  return _specIntegratorConfig ;
}


////////////////////////////////////////////////////////////////////////////////
/// Returns the specialized integrator configuration for _this_ RooAbsReal.
/// If this object has no specialized configuration, a null pointer is returned,
/// unless createOnTheFly is kTRUE in which case a clone of the default integrator
/// configuration is created, installed as specialized configuration, and returned

RooNumIntConfig* RooAbsReal::specialIntegratorConfig(Bool_t createOnTheFly)
{
  if (!_specIntegratorConfig && createOnTheFly) {
    _specIntegratorConfig = new RooNumIntConfig(*defaultIntegratorConfig()) ;
  }
  return _specIntegratorConfig ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return the numeric integration configuration used for this object. If
/// a specialized configuration was associated with this object, that configuration
/// is returned, otherwise the default configuration for all RooAbsReals is returned

const RooNumIntConfig* RooAbsReal::getIntegratorConfig() const
{
  const RooNumIntConfig* config = specialIntegratorConfig() ;
  if (config) return config ;
  return defaultIntegratorConfig() ;
}


////////////////////////////////////////////////////////////////////////////////
/// Return the numeric integration configuration used for this object. If
/// a specialized configuration was associated with this object, that configuration
/// is returned, otherwise the default configuration for all RooAbsReals is returned

RooNumIntConfig* RooAbsReal::getIntegratorConfig()
{
  RooNumIntConfig* config = specialIntegratorConfig() ;
  if (config) return config ;
  return defaultIntegratorConfig() ;
}



////////////////////////////////////////////////////////////////////////////////
/// Set the given integrator configuration as default numeric integration
/// configuration for this object

void RooAbsReal::setIntegratorConfig(const RooNumIntConfig& config)
{
  if (_specIntegratorConfig) {
    delete _specIntegratorConfig ;
  }
  _specIntegratorConfig = new RooNumIntConfig(config) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Remove the specialized numeric integration configuration associated
/// with this object

void RooAbsReal::setIntegratorConfig()
{
  if (_specIntegratorConfig) {
    delete _specIntegratorConfig ;
  }
  _specIntegratorConfig = 0 ;
}




////////////////////////////////////////////////////////////////////////////////
/// Interface function to force use of a given set of observables
/// to interpret function value. Needed for functions or p.d.f.s
/// whose shape depends on the choice of normalization such as
/// RooAddPdf

void RooAbsReal::selectNormalization(const RooArgSet*, Bool_t)
{
}




////////////////////////////////////////////////////////////////////////////////
/// Interface function to force use of a given normalization range
/// to interpret function value. Needed for functions or p.d.f.s
/// whose shape depends on the choice of normalization such as
/// RooAddPdf

void RooAbsReal::selectNormalizationRange(const char*, Bool_t)
{
}



////////////////////////////////////////////////////////////////////////////////
/// Activate cache validation mode

void RooAbsReal::setCacheCheck(Bool_t flag)
{
  _cacheCheck = flag ;
}



////////////////////////////////////////////////////////////////////////////////
/// Advertise capability to determine maximum value of function for given set of
/// observables. If no direct generator method is provided, this information
/// will assist the accept/reject generator to operate more efficiently as
/// it can skip the initial trial sampling phase to empirically find the function
/// maximum

Int_t RooAbsReal::getMaxVal(const RooArgSet& /*vars*/) const
{
  return 0 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return maximum value for set of observables identified by code assigned
/// in getMaxVal

Double_t RooAbsReal::maxVal(Int_t /*code*/) const
{
  assert(1) ;
  return 0 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Interface to insert remote error logging messages received by RooRealMPFE into current error loggin stream

void RooAbsReal::logEvalError(const RooAbsReal* originator, const char* origName, const char* message, const char* serverValueString)
{
  if (_evalErrorMode==Ignore) {
    return ;
  }

  if (_evalErrorMode==CountErrors) {
    _evalErrorCount++ ;
    return ;
  }

  static Bool_t inLogEvalError = kFALSE ;

  if (inLogEvalError) {
    return ;
  }
  inLogEvalError = kTRUE ;

  EvalError ee ;
  ee.setMessage(message) ;

  if (serverValueString) {
    ee.setServerValues(serverValueString) ;
  }

  if (_evalErrorMode==PrintErrors) {
   oocoutE((TObject*)0,Eval) << "RooAbsReal::logEvalError(" << "<STATIC>" << ") evaluation error, " << endl
         << " origin       : " << origName << endl
         << " message      : " << ee._msg << endl
         << " server values: " << ee._srvval << endl ;
  } else if (_evalErrorMode==CollectErrors) {
    _evalErrorList[originator].first = origName ;
    _evalErrorList[originator].second.push_back(ee) ;
  }


  inLogEvalError = kFALSE ;
}



////////////////////////////////////////////////////////////////////////////////
/// Log evaluation error message. Evaluation errors may be routed through a different
/// protocol than generic RooFit warning message (which go straight through RooMsgService)
/// because evaluation errors can occur in very large numbers in the use of likelihood
/// evaluations. In logEvalError mode, controlled by global method enableEvalErrorLogging()
/// messages reported through this function are not printed but all stored in a list,
/// along with server values at the time of reporting. Error messages logged in this
/// way can be printed in a structured way, eliminating duplicates and with the ability
/// to truncate the list by printEvalErrors. This is the standard mode of error logging
/// during MINUIT operations. If enableEvalErrorLogging() is false, all errors
/// reported through this method are passed for immediate printing through RooMsgService.
/// A string with server names and values is constructed automatically for error logging
/// purposes, unless a custom string with similar information is passed as argument.

void RooAbsReal::logEvalError(const char* message, const char* serverValueString) const
{
  if (_evalErrorMode==Ignore) {
    return ;
  }

  if (_evalErrorMode==CountErrors) {
    _evalErrorCount++ ;
    return ;
  }

  static Bool_t inLogEvalError = kFALSE ;

  if (inLogEvalError) {
    return ;
  }
  inLogEvalError = kTRUE ;

  EvalError ee ;
  ee.setMessage(message) ;

  if (serverValueString) {
    ee.setServerValues(serverValueString) ;
  } else {
    string srvval ;
    ostringstream oss ;
    Bool_t first(kTRUE) ;
    for (Int_t i=0 ; i<numProxies() ; i++) {
      RooAbsProxy* p = getProxy(i) ;
      if (!p) continue ;
      //if (p->name()[0]=='!') continue ;
      if (first) {
   first=kFALSE ;
      } else {
   oss << ", " ;
      }
      p->print(oss,kTRUE) ;
    }
    ee.setServerValues(oss.str().c_str()) ;
  }

  ostringstream oss2 ;
  printStream(oss2,kName|kClassName|kArgs,kInline)  ;

  if (_evalErrorMode==PrintErrors) {
   coutE(Eval) << "RooAbsReal::logEvalError(" << GetName() << ") evaluation error, " << endl
          << " origin       : " << oss2.str() << endl
          << " message      : " << ee._msg << endl
          << " server values: " << ee._srvval << endl ;
  } else if (_evalErrorMode==CollectErrors) {
    if (_evalErrorList[this].second.size() >= 2048) {
       // avoid overflowing the error list, so if there are very many, print
       // the oldest one first, and pop it off the list
       const EvalError& oee = _evalErrorList[this].second.front();
       // print to debug stream, since these would normally be suppressed, and
       // we do not want to increase the error count in the message service...
       ccoutD(Eval) << "RooAbsReal::logEvalError(" << GetName()
              << ") delayed evaluation error, " << endl
                   << " origin       : " << oss2.str() << endl
                   << " message      : " << oee._msg << endl
                   << " server values: " << oee._srvval << endl ;
       _evalErrorList[this].second.pop_front();
    }
    _evalErrorList[this].first = oss2.str().c_str() ;
    _evalErrorList[this].second.push_back(ee) ;
  }

  inLogEvalError = kFALSE ;
  //coutE(Tracing) << "RooAbsReal::logEvalError(" << GetName() << ") message = " << message << endl ;
}




////////////////////////////////////////////////////////////////////////////////
/// Clear the stack of evaluation error messages

void RooAbsReal::clearEvalErrorLog()
{
  if (_evalErrorMode==PrintErrors) {
    return ;
  } else if (_evalErrorMode==CollectErrors) {
    _evalErrorList.clear() ;
  } else {
    _evalErrorCount = 0 ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Print all outstanding logged evaluation error on the given ostream. If maxPerNode
/// is zero, only the number of errors for each source (object with unique name) is listed.
/// If maxPerNode is greater than zero, up to maxPerNode detailed error messages are shown
/// per source of errors. A truncation message is shown if there were more errors logged
/// than shown.

void RooAbsReal::printEvalErrors(ostream& os, Int_t maxPerNode)
{
  if (_evalErrorMode == CountErrors) {
    os << _evalErrorCount << " errors counted" << endl ;
  }

  if (maxPerNode<0) return ;

  map<const RooAbsArg*,pair<string,list<EvalError> > >::iterator iter = _evalErrorList.begin() ;

  for(;iter!=_evalErrorList.end() ; ++iter) {
    if (maxPerNode==0) {

      // Only print node name with total number of errors
      os << iter->second.first ;
      //iter->first->printStream(os,kName|kClassName|kArgs,kInline)  ;
      os << " has " << iter->second.second.size() << " errors" << endl ;

    } else {

      // Print node name and details of 'maxPerNode' errors
      os << iter->second.first << endl ;
      //iter->first->printStream(os,kName|kClassName|kArgs,kSingleLine) ;

      Int_t i(0) ;
      std::list<EvalError>::iterator iter2 = iter->second.second.begin() ;
      for(;iter2!=iter->second.second.end() ; ++iter2, i++) {
   os << "     " << iter2->_msg << " @ " << iter2->_srvval << endl ;
   if (i>maxPerNode) {
     os << "    ... (remaining " << iter->second.second.size() - maxPerNode << " messages suppressed)" << endl ;
     break ;
   }
      }
    }
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Return the number of logged evaluation errors since the last clearing.

Int_t RooAbsReal::numEvalErrors()
{
  if (_evalErrorMode==CountErrors) {
    return _evalErrorCount ;
  }

  Int_t ntot(0) ;
  map<const RooAbsArg*,pair<string,list<EvalError> > >::iterator iter = _evalErrorList.begin() ;
  for(;iter!=_evalErrorList.end() ; ++iter) {
    ntot += iter->second.second.size() ;
  }
  return ntot ;
}



////////////////////////////////////////////////////////////////////////////////
/// Fix the interpretation of the coefficient of any RooAddPdf component in
/// the expression tree headed by this object to the given set of observables.
///
/// If the force flag is false, the normalization choice is only fixed for those
/// RooAddPdf components that have the default 'automatic' interpretation of
/// coefficients (i.e. the interpretation is defined by the observables passed
/// to getVal()). If force is true, also RooAddPdf that already have a fixed
/// interpretation are changed to a new fixed interpretation.

void RooAbsReal::fixAddCoefNormalization(const RooArgSet& addNormSet, Bool_t force)
{
  RooArgSet* compSet = getComponents() ;
  TIterator* iter = compSet->createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
    if (pdf) {
      if (addNormSet.getSize()>0) {
   pdf->selectNormalization(&addNormSet,force) ;
      } else {
   pdf->selectNormalization(0,force) ;
      }
    }
  }
  delete iter ;
  delete compSet ;
}



////////////////////////////////////////////////////////////////////////////////
/// Fix the interpretation of the coefficient of any RooAddPdf component in
/// the expression tree headed by this object to the given set of observables.
///
/// If the force flag is false, the normalization range choice is only fixed for those
/// RooAddPdf components that currently use the default full domain to interpret their
/// coefficients. If force is true, also RooAddPdf that already have a fixed
/// interpretation range are changed to a new fixed interpretation range.

void RooAbsReal::fixAddCoefRange(const char* rangeName, Bool_t force)
{
  RooArgSet* compSet = getComponents() ;
  TIterator* iter = compSet->createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
    if (pdf) {
      pdf->selectNormalizationRange(rangeName,force) ;
    }
  }
  delete iter ;
  delete compSet ;
}



////////////////////////////////////////////////////////////////////////////////
/// Interface method for function objects to indicate their prefferred order of observables
/// for scanning their values into a (multi-dimensional) histogram or RooDataSet. The observables
/// to be ordered are offered in argument 'obs' and should be copied in their preferred
/// order into argument 'orderdObs', This default implementation indicates no preference
/// and copies the original order of 'obs' into 'orderedObs'

void RooAbsReal::preferredObservableScanOrder(const RooArgSet& obs, RooArgSet& orderedObs) const
{
  // Dummy implementation, do nothing
  orderedObs.removeAll() ;
  orderedObs.add(obs) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Create a running integral over this function, i.e. given a f(x), create an object
/// representing 'int[x_lo,x] f(x_prime) dx_prime'

RooAbsReal* RooAbsReal::createRunningIntegral(const RooArgSet& iset, const RooArgSet& nset)
{
  return createRunningIntegral(iset,RooFit::SupNormSet(nset)) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Create an object that represents the running integral of the function over one or more observables listed in iset, i.e.
///
///   int[x_lo,x] f(x_prime) dx_prime
///
/// The actual integration calculation is only performed when the return object is evaluated. The name
/// of the integral object is automatically constructed from the name of the input function, the variables
/// it integrates and the range integrates over. The default strategy to calculate the running integrals is
///
///   - If the integrand (this object) supports analytical integration, construct an integral object
///     that calculate the running integrals value by calculating the analytical integral each
///     time the running integral object is evaluated
///
///   - If the integrand (this object) requires numeric integration to construct the running integral
///     create an object of class RooNumRunningInt which first samples the entire function and integrates
///     the sampled function numerically. This method has superior performance as there is no need to
///     perform a full (numeric) integration for each evaluation of the running integral object, but
///     only when one of its parameters has changed.
///
/// The choice of strategy can be changed with the ScanAll() argument, which forces the use of the
/// scanning technique implemented in RooNumRunningInt for all use cases, and with the ScanNone()
/// argument which forces the 'integrate each evaluation' technique for all use cases. The sampling
/// granularity for the scanning technique can be controlled with the ScanParameters technique
/// which allows to specify the number of samples to be taken, and to which order the resulting
/// running integral should be interpolated. The default values are 1000 samples and 2nd order
/// interpolation.
///
/// The following named arguments are accepted
///
/// SupNormSet(const RooArgSet&)         -- Observables over which should be normalized _in_addition_ to the
///                                         integration observables
/// ScanParameters(Int_t nbins,          -- Parameters for scanning technique of making CDF: number
///                Int_t intOrder)          of sampled bins and order of interpolation applied on numeric cdf
/// ScanNum()                            -- Apply scanning technique if cdf integral involves numeric integration
/// ScanAll()                            -- Always apply scanning technique
/// ScanNone()                           -- Never apply scanning technique

RooAbsReal* RooAbsReal::createRunningIntegral(const RooArgSet& iset, const RooCmdArg& arg1, const RooCmdArg& arg2,
             const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
             const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
{
  // Define configuration for this method
  RooCmdConfig pc(Form("RooAbsReal::createRunningIntegral(%s)",GetName())) ;
  pc.defineObject("supNormSet","SupNormSet",0,0) ;
  pc.defineInt("numScanBins","ScanParameters",0,1000) ;
  pc.defineInt("intOrder","ScanParameters",1,2) ;
  pc.defineInt("doScanNum","ScanNum",0,1) ;
  pc.defineInt("doScanAll","ScanAll",0,0) ;
  pc.defineInt("doScanNon","ScanNone",0,0) ;
  pc.defineMutex("ScanNum","ScanAll","ScanNone") ;

  // Process & check varargs
  pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
  if (!pc.ok(kTRUE)) {
    return 0 ;
  }

  // Extract values from named arguments
  const RooArgSet* snset = static_cast<const RooArgSet*>(pc.getObject("supNormSet",0)) ;
  RooArgSet nset ;
  if (snset) {
    nset.add(*snset) ;
  }
  Int_t numScanBins = pc.getInt("numScanBins") ;
  Int_t intOrder = pc.getInt("intOrder") ;
  Int_t doScanNum = pc.getInt("doScanNum") ;
  Int_t doScanAll = pc.getInt("doScanAll") ;
  Int_t doScanNon = pc.getInt("doScanNon") ;

  // If scanning technique is not requested make integral-based cdf and return
  if (doScanNon) {
    return createIntRI(iset,nset) ;
  }
  if (doScanAll) {
    return createScanRI(iset,nset,numScanBins,intOrder) ;
  }
  if (doScanNum) {
    RooRealIntegral* tmp = (RooRealIntegral*) createIntegral(iset) ;
    Int_t isNum= (tmp->numIntRealVars().getSize()==1) ;
    delete tmp ;

    if (isNum) {
      coutI(NumIntegration) << "RooAbsPdf::createRunningIntegral(" << GetName() << ") integration over observable(s) " << iset << " involves numeric integration," << endl
             << "      constructing cdf though numeric integration of sampled pdf in " << numScanBins << " bins and applying order "
             << intOrder << " interpolation on integrated histogram." << endl
             << "      To override this choice of technique use argument ScanNone(), to change scan parameters use ScanParameters(nbins,order) argument" << endl ;
    }

    return isNum ? createScanRI(iset,nset,numScanBins,intOrder) : createIntRI(iset,nset) ;
  }
  return 0 ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for createRunningIntegral that construct an object
/// implementing the numeric scanning technique for calculating the running integral

RooAbsReal* RooAbsReal::createScanRI(const RooArgSet& iset, const RooArgSet& nset, Int_t numScanBins, Int_t intOrder)
{
  string name = string(GetName()) + "_NUMRUNINT_" + integralNameSuffix(iset,&nset).Data() ;
  RooRealVar* ivar = (RooRealVar*) iset.first() ;
  ivar->setBins(numScanBins,"numcdf") ;
  RooNumRunningInt* ret = new RooNumRunningInt(name.c_str(),name.c_str(),*this,*ivar,"numrunint") ;
  ret->setInterpolationOrder(intOrder) ;
  return ret ;
}



////////////////////////////////////////////////////////////////////////////////
/// Utility function for createRunningIntegral that construct an
/// object implementing the standard (analytical) integration
/// technique for calculating the running integral

RooAbsReal* RooAbsReal::createIntRI(const RooArgSet& iset, const RooArgSet& nset)
{
  // Make list of input arguments keeping only RooRealVars
  RooArgList ilist ;
  TIterator* iter2 = iset.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter2->Next())) {
    if (dynamic_cast<RooRealVar*>(arg)) {
      ilist.add(*arg) ;
    } else {
      coutW(InputArguments) << "RooAbsPdf::createRunningIntegral(" << GetName() << ") WARNING ignoring non-RooRealVar input argument " << arg->GetName() << endl ;
    }
  }
  delete iter2 ;

  RooArgList cloneList ;
  RooArgList loList ;
  RooArgSet clonedBranchNodes ;

  // Setup customizer that stores all cloned branches in our non-owning list
  RooCustomizer cust(*this,"cdf") ;
  cust.setCloneBranchSet(clonedBranchNodes) ;
  cust.setOwning(kFALSE) ;

  // Make integration observable x_prime for each observable x as well as an x_lowbound
  TIterator* iter = ilist.createIterator() ;
  RooRealVar* rrv ;
  while((rrv=(RooRealVar*)iter->Next())) {

    // Make clone x_prime of each c.d.f observable x represening running integral
    RooRealVar* cloneArg = (RooRealVar*) rrv->clone(Form("%s_prime",rrv->GetName())) ;
    cloneList.add(*cloneArg) ;
    cust.replaceArg(*rrv,*cloneArg) ;

    // Make clone x_lowbound of each c.d.f observable representing low bound of x
    RooRealVar* cloneLo = (RooRealVar*) rrv->clone(Form("%s_lowbound",rrv->GetName())) ;
    cloneLo->setVal(rrv->getMin()) ;
    loList.add(*cloneLo) ;

    // Make parameterized binning from [x_lowbound,x] for each x_prime
    RooParamBinning pb(*cloneLo,*rrv,100) ;
    cloneArg->setBinning(pb,"CDF") ;

  }
  delete iter ;

  RooAbsReal* tmp = (RooAbsReal*) cust.build() ;

  // Construct final normalization set for c.d.f = integrated observables + any extra specified by user
  RooArgSet finalNset(nset) ;
  finalNset.add(cloneList,kTRUE) ;
  RooAbsReal* cdf = tmp->createIntegral(cloneList,finalNset,"CDF") ;

  // Transfer ownership of cloned items to top-level c.d.f object
  cdf->addOwnedComponents(*tmp) ;
  cdf->addOwnedComponents(cloneList) ;
  cdf->addOwnedComponents(loList) ;

  return cdf ;
}


////////////////////////////////////////////////////////////////////////////////
/// Return a RooFunctor object bound to this RooAbsReal with given definition of observables
/// and parameters

RooFunctor* RooAbsReal::functor(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
{
  RooArgSet* realObs = getObservables(obs) ;
  if (realObs->getSize() != obs.getSize()) {
    coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified observables are not variables of this p.d.f" << endl ;
    delete realObs ;
    return 0 ;
  }
  RooArgSet* realPars = getObservables(pars) ;
  if (realPars->getSize() != pars.getSize()) {
    coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified parameters are not variables of this p.d.f" << endl ;
    delete realPars ;
    return 0 ;
  }
  delete realObs ;
  delete realPars ;

  return new RooFunctor(*this,obs,pars,nset) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return a ROOT TF1,2,3 object bound to this RooAbsReal with given definition of observables
/// and parameters

TF1* RooAbsReal::asTF(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
{
  // Check that specified input are indeed variables of this function
  RooArgSet* realObs = getObservables(obs) ;
  if (realObs->getSize() != obs.getSize()) {
    coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified observables are not variables of this p.d.f" << endl ;
    delete realObs ;
    return 0 ;
  }
  RooArgSet* realPars = getObservables(pars) ;
  if (realPars->getSize() != pars.getSize()) {
    coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified parameters are not variables of this p.d.f" << endl ;
    delete realPars ;
    return 0 ;
  }
  delete realObs ;
  delete realPars ;

  // Check that all obs and par are of type RooRealVar
  for (int i=0 ; i<obs.getSize() ; i++) {
    if (dynamic_cast<RooRealVar*>(obs.at(i))==0) {
      coutE(ObjectHandling) << "RooAbsReal::asTF(" << GetName() << ") ERROR: proposed observable " << obs.at(0)->GetName() << " is not of type RooRealVar" << endl ;
      return 0 ;
    }
  }
  for (int i=0 ; i<pars.getSize() ; i++) {
    if (dynamic_cast<RooRealVar*>(pars.at(i))==0) {
      coutE(ObjectHandling) << "RooAbsReal::asTF(" << GetName() << ") ERROR: proposed parameter " << pars.at(0)->GetName() << " is not of type RooRealVar" << endl ;
      return 0 ;
    }
  }

  // Create functor and TFx of matching dimension
  TF1* tf=0 ;
  RooFunctor* f ;
  switch(obs.getSize()) {
  case 1: {
    RooRealVar* x = (RooRealVar*)obs.at(0) ;
    f = functor(obs,pars,nset) ;
    tf = new TF1(GetName(),f,x->getMin(),x->getMax(),pars.getSize()) ;
    break ;
  }
  case 2: {
    RooRealVar* x = (RooRealVar*)obs.at(0) ;
    RooRealVar* y = (RooRealVar*)obs.at(1) ;
    f = functor(obs,pars,nset) ;
    tf = new TF2(GetName(),f,x->getMin(),x->getMax(),y->getMin(),y->getMax(),pars.getSize()) ;
    break ;
  }
  case 3: {
    RooRealVar* x = (RooRealVar*)obs.at(0) ;
    RooRealVar* y = (RooRealVar*)obs.at(1) ;
    RooRealVar* z = (RooRealVar*)obs.at(2) ;
    f = functor(obs,pars,nset) ;
    tf = new TF3(GetName(),f,x->getMin(),x->getMax(),y->getMin(),y->getMax(),z->getMin(),z->getMax(),pars.getSize()) ;
    break ;
  }
  default:
    coutE(InputArguments) << "RooAbsReal::asTF(" << GetName() << ") ERROR: " << obs.getSize()
           << " observables specified, but a ROOT TFx can only have  1,2 or 3 observables" << endl ;
    return 0 ;
  }

  // Set initial parameter values of TFx to those of RooRealVars
  for (int i=0 ; i<pars.getSize() ; i++) {
    RooRealVar* p = (RooRealVar*) pars.at(i) ;
    tf->SetParameter(i,p->getVal()) ;
    tf->SetParName(i,p->GetName()) ;
    //tf->SetParLimits(i,p->getMin(),p->getMax()) ;
  }

  return tf ;
}


////////////////////////////////////////////////////////////////////////////////
/// Return function representing first, second or third order derivative of this function

RooDerivative* RooAbsReal::derivative(RooRealVar& obs, Int_t order, Double_t eps)
{
  string name=Form("%s_DERIV_%s",GetName(),obs.GetName()) ;
  string title=Form("Derivative of %s w.r.t %s ",GetName(),obs.GetName()) ;
  return new RooDerivative(name.c_str(),title.c_str(),*this,obs,order,eps) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return function representing first, second or third order derivative of this function

RooDerivative* RooAbsReal::derivative(RooRealVar& obs, const RooArgSet& normSet, Int_t order, Double_t eps)
{
  string name=Form("%s_DERIV_%s",GetName(),obs.GetName()) ;
  string title=Form("Derivative of %s w.r.t %s ",GetName(),obs.GetName()) ;
  return new RooDerivative(name.c_str(),title.c_str(),*this,obs,normSet,order,eps) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return function representing moment of function of given order. If central is
/// true, the central moment is given <(x-<x>)^2>

RooAbsMoment* RooAbsReal::moment(RooRealVar& obs, Int_t order, Bool_t central, Bool_t takeRoot)
{
  string name=Form("%s_MOMENT_%d%s_%s",GetName(),order,(central?"C":""),obs.GetName()) ;
  string title=Form("%sMoment of order %d of %s w.r.t %s ",(central?"Central ":""),order,GetName(),obs.GetName()) ;
  if (order==1) return new RooFirstMoment(name.c_str(),title.c_str(),*this,obs) ;
  if (order==2) return new RooSecondMoment(name.c_str(),title.c_str(),*this,obs,central,takeRoot) ;
  return new RooMoment(name.c_str(),title.c_str(),*this,obs,order,central,takeRoot) ;
}


////////////////////////////////////////////////////////////////////////////////
/// Return function representing moment of p.d.f (normalized w.r.t given observables) of given order. If central is
/// true, the central moment is given <(x-<x>)^2>. If intNormObs is true, the moment of the function integrated over
/// all normalization observables is returned.

RooAbsMoment* RooAbsReal::moment(RooRealVar& obs, const RooArgSet& normObs, Int_t order, Bool_t central, Bool_t takeRoot, Bool_t intNormObs)
{
  string name=Form("%s_MOMENT_%d%s_%s",GetName(),order,(central?"C":""),obs.GetName()) ;
  string title=Form("%sMoment of order %d of %s w.r.t %s ",(central?"Central ":""),order,GetName(),obs.GetName()) ;

  if (order==1) return new RooFirstMoment(name.c_str(),title.c_str(),*this,obs,normObs,intNormObs) ;
  if (order==2) return new RooSecondMoment(name.c_str(),title.c_str(),*this,obs,normObs,central,takeRoot,intNormObs) ;
  return new RooMoment(name.c_str(),title.c_str(),*this,obs,normObs,order,central,takeRoot,intNormObs) ;
}



////////////////////////////////////////////////////////////////////////////////
///
/// Return value of x (in range xmin,xmax) at which function equals yval.
/// (Calculation is performed with Brent root finding algorithm)

Double_t RooAbsReal::findRoot(RooRealVar& x, Double_t xmin, Double_t xmax, Double_t yval)
{
  Double_t result(0) ;
  RooBrentRootFinder(RooRealBinding(*this,x)).findRoot(result,xmin,xmax,yval) ;
  return result ;
}




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

RooGenFunction* RooAbsReal::iGenFunction(RooRealVar& x, const RooArgSet& nset)
{
  return new RooGenFunction(*this,x,RooArgList(),nset.getSize()>0?nset:RooArgSet(x)) ;
}



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

RooMultiGenFunction* RooAbsReal::iGenFunction(const RooArgSet& observables, const RooArgSet& nset)
{
  return new RooMultiGenFunction(*this,observables,RooArgList(),nset.getSize()>0?nset:observables) ;
}




////////////////////////////////////////////////////////////////////////////////
/// Perform a chi^2 fit to given histogram By default the fit is executed through the MINUIT
/// commands MIGRAD, HESSE in succession
///
/// The following named arguments are supported
///
/// Options to control construction of -log(L)
/// ------------------------------------------
/// Range(const char* name)         -- Fit only data inside range with given name
/// Range(Double_t lo, Double_t hi) -- Fit only data inside given range. A range named "fit" is created on the fly on all observables.
///                                    Multiple comma separated range names can be specified.
/// NumCPU(int num)                 -- Parallelize NLL calculation on num CPUs
/// Optimize(Bool_t flag)           -- Activate constant term optimization (on by default)
///
/// Options to control flow of fit procedure
/// ----------------------------------------
/// InitialHesse(Bool_t flag)      -- Flag controls if HESSE before MIGRAD as well, off by default
/// Hesse(Bool_t flag)             -- Flag controls if HESSE is run after MIGRAD, on by default
/// Minos(Bool_t flag)             -- Flag controls if MINOS is run after HESSE, on by default
/// Minos(const RooArgSet& set)    -- Only run MINOS on given subset of arguments
/// Save(Bool_t flag)              -- Flac controls if RooFitResult object is produced and returned, off by default
/// Strategy(Int_t flag)           -- Set Minuit strategy (0 through 2, default is 1)
/// FitOptions(const char* optStr) -- Steer fit with classic options string (for backward compatibility). Use of this option
///                                   excludes use of any of the new style steering options.
///
/// Options to control informational output
/// ---------------------------------------
/// Verbose(Bool_t flag)           -- Flag controls if verbose output is printed (NLL, parameter changes during fit
/// Timer(Bool_t flag)             -- Time CPU and wall clock consumption of fit steps, off by default
/// PrintLevel(Int_t level)        -- Set Minuit print level (-1 through 3, default is 1). At -1 all RooFit informational
///                                   messages are suppressed as well
/// Warnings(Bool_t flag)          -- Enable or disable MINUIT warnings (enabled by default)
/// PrintEvalErrors(Int_t numErr)  -- Control number of p.d.f evaluation errors printed per likelihood evaluation. A negative
///                                   value suppress output completely, a zero value will only print the error count per p.d.f component,
///                                   a positive value is will print details of each error up to numErr messages per p.d.f component.
///
///

RooFitResult* RooAbsReal::chi2FitTo(RooDataHist& data, const RooCmdArg& arg1,  const RooCmdArg& arg2,
                const RooCmdArg& arg3,  const RooCmdArg& arg4, const RooCmdArg& arg5,
                const RooCmdArg& arg6,  const RooCmdArg& arg7, const RooCmdArg& arg8)
{
  RooLinkedList l ;
  l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
  l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
  l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
  l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;
  return chi2FitTo(data,l) ;

}



////////////////////////////////////////////////////////////////////////////////
/// Internal back-end function to steer chi2 fits

RooFitResult* RooAbsReal::chi2FitTo(RooDataHist& data, const RooLinkedList& cmdList)
{
  // Select the pdf-specific commands
  RooCmdConfig pc(Form("RooAbsPdf::chi2FitTo(%s)",GetName())) ;

  // Pull arguments to be passed to chi2 construction from list
  RooLinkedList fitCmdList(cmdList) ;
  RooLinkedList chi2CmdList = pc.filterCmdList(fitCmdList,"Range,RangeWithName,NumCPU,Optimize") ;

  RooAbsReal* chi2 = createChi2(data,chi2CmdList) ;
  RooFitResult* ret = chi2FitDriver(*chi2,fitCmdList) ;

  // Cleanup
  delete chi2 ;
  return ret ;
}




////////////////////////////////////////////////////////////////////////////////
/// Create a chi-2 from a histogram and this function.
///
/// The following named arguments are supported
///
///  Options to control construction of the chi^2
///  ------------------------------------------
///  DataError(RooAbsData::ErrorType)  -- Choose between Poisson errors and Sum-of-weights errors
///  NumCPU(Int_t)                     -- Activate parallel processing feature on N processes
///  Range()                           -- Calculate Chi2 only in selected region

RooAbsReal* RooAbsReal::createChi2(RooDataHist& data, const RooCmdArg& arg1,  const RooCmdArg& arg2,
               const RooCmdArg& arg3,  const RooCmdArg& arg4, const RooCmdArg& arg5,
               const RooCmdArg& arg6,  const RooCmdArg& arg7, const RooCmdArg& arg8)
{
  string name = Form("chi2_%s_%s",GetName(),data.GetName()) ;

  return new RooChi2Var(name.c_str(),name.c_str(),*this,data,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
}




////////////////////////////////////////////////////////////////////////////////
/// Internal back-end function to create a chi2

RooAbsReal* RooAbsReal::createChi2(RooDataHist& data, const RooLinkedList& cmdList)
{
  // Fill array of commands
  const RooCmdArg* cmds[8] ;
  TIterator* iter = cmdList.MakeIterator() ;
  Int_t i(0) ;
  RooCmdArg* arg ;
  while((arg=(RooCmdArg*)iter->Next())) {
    cmds[i++] = arg ;
  }
  for (;i<8 ; i++) {
    cmds[i] = &RooCmdArg::none() ;
  }
  delete iter ;

  return createChi2(data,*cmds[0],*cmds[1],*cmds[2],*cmds[3],*cmds[4],*cmds[5],*cmds[6],*cmds[7]) ;

}





////////////////////////////////////////////////////////////////////////////////
/// Create a chi-2 from a series of x and y value stored in a dataset.
/// The y values can either be the event weights, or can be another column designated
/// by the YVar() argument. The y value must have errors defined for the chi-2 to
/// be well defined.
///
/// The following named arguments are supported
///
/// Options to control construction of the chi^2
/// ------------------------------------------
/// YVar(RooRealVar& yvar)          -- Designate given column in dataset as Y value
/// Integrate(Bool_t flag)          -- Integrate function over range specified by X errors
///                                    rather than take value at bin center.
///
/// Options to control flow of fit procedure
/// ----------------------------------------
/// InitialHesse(Bool_t flag)      -- Flag controls if HESSE before MIGRAD as well, off by default
/// Hesse(Bool_t flag)             -- Flag controls if HESSE is run after MIGRAD, on by default
/// Minos(Bool_t flag)             -- Flag controls if MINOS is run after HESSE, on by default
/// Minos(const RooArgSet& set)    -- Only run MINOS on given subset of arguments
/// Save(Bool_t flag)              -- Flac controls if RooFitResult object is produced and returned, off by default
/// Strategy(Int_t flag)           -- Set Minuit strategy (0 through 2, default is 1)
/// FitOptions(const char* optStr) -- Steer fit with classic options string (for backward compatibility). Use of this option
///                                   excludes use of any of the new style steering options.
///
/// Options to control informational output
/// ---------------------------------------
/// Verbose(Bool_t flag)           -- Flag controls if verbose output is printed (NLL, parameter changes during fit
/// Timer(Bool_t flag)             -- Time CPU and wall clock consumption of fit steps, off by default
/// PrintLevel(Int_t level)        -- Set Minuit print level (-1 through 3, default is 1). At -1 all RooFit informational
///                                   messages are suppressed as well
/// Warnings(Bool_t flag)          -- Enable or disable MINUIT warnings (enabled by default)
/// PrintEvalErrors(Int_t numErr)  -- Control number of p.d.f evaluation errors printed per likelihood evaluation. A negative
///                                   value suppress output completely, a zero value will only print the error count per p.d.f component,
///                                   a positive value is will print details of each error up to numErr messages per p.d.f component.

RooFitResult* RooAbsReal::chi2FitTo(RooDataSet& xydata, const RooCmdArg& arg1,  const RooCmdArg& arg2,
                  const RooCmdArg& arg3,  const RooCmdArg& arg4, const RooCmdArg& arg5,
                  const RooCmdArg& arg6,  const RooCmdArg& arg7, const RooCmdArg& arg8)
{
  RooLinkedList l ;
  l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
  l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
  l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
  l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;
  return chi2FitTo(xydata,l) ;
}




////////////////////////////////////////////////////////////////////////////////
/// Internal back-end function to steer chi2 fits

RooFitResult* RooAbsReal::chi2FitTo(RooDataSet& xydata, const RooLinkedList& cmdList)
{
  // Select the pdf-specific commands
  RooCmdConfig pc(Form("RooAbsPdf::chi2FitTo(%s)",GetName())) ;

  // Pull arguments to be passed to chi2 construction from list
  RooLinkedList fitCmdList(cmdList) ;
  RooLinkedList chi2CmdList = pc.filterCmdList(fitCmdList,"YVar,Integrate") ;

  RooAbsReal* xychi2 = createChi2(xydata,chi2CmdList) ;
  RooFitResult* ret = chi2FitDriver(*xychi2,fitCmdList) ;

  // Cleanup
  delete xychi2 ;
  return ret ;
}




////////////////////////////////////////////////////////////////////////////////
/// Create a chi-2 from a series of x and y value stored in a dataset.
/// The y values can either be the event weights (default), or can be another column designated
/// by the YVar() argument. The y value must have errors defined for the chi-2 to
/// be well defined.
///
/// The following named arguments are supported
///
/// Options to control construction of the chi^2
/// ------------------------------------------
/// YVar(RooRealVar& yvar)          -- Designate given column in dataset as Y value
/// Integrate(Bool_t flag)          -- Integrate function over range specified by X errors
///                                    rather than take value at bin center.
///

RooAbsReal* RooAbsReal::createChi2(RooDataSet& data, const RooCmdArg& arg1,  const RooCmdArg& arg2,
                 const RooCmdArg& arg3,  const RooCmdArg& arg4, const RooCmdArg& arg5,
                 const RooCmdArg& arg6,  const RooCmdArg& arg7, const RooCmdArg& arg8)
{
  RooLinkedList l ;
  l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
  l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
  l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
  l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;
  return createChi2(data,l) ;
}



////////////////////////////////////////////////////////////////////////////////
/// Internal back-end function to create a chi^2 from a function and a dataset

RooAbsReal* RooAbsReal::createChi2(RooDataSet& data, const RooLinkedList& cmdList)
{
  // Select the pdf-specific commands
  RooCmdConfig pc(Form("RooAbsPdf::fitTo(%s)",GetName())) ;

  pc.defineInt("integrate","Integrate",0,0) ;
  pc.defineObject("yvar","YVar",0,0) ;

  // Process and check varargs
  pc.process(cmdList) ;
  if (!pc.ok(kTRUE)) {
    return 0 ;
  }

  // Decode command line arguments
  Bool_t integrate = pc.getInt("integrate") ;
  RooRealVar* yvar = (RooRealVar*) pc.getObject("yvar") ;

  string name = Form("chi2_%s_%s",GetName(),data.GetName()) ;

  if (yvar) {
    return new RooXYChi2Var(name.c_str(),name.c_str(),*this,data,*yvar,integrate) ;
  } else {
    return new RooXYChi2Var(name.c_str(),name.c_str(),*this,data,integrate) ;
  }
}






////////////////////////////////////////////////////////////////////////////////
/// Internal driver function for chi2 fits

RooFitResult* RooAbsReal::chi2FitDriver(RooAbsReal& fcn, RooLinkedList& cmdList)
{
  // Select the pdf-specific commands
  RooCmdConfig pc(Form("RooAbsPdf::chi2FitDriver(%s)",GetName())) ;

  pc.defineString("fitOpt","FitOptions",0,"") ;

  pc.defineInt("optConst","Optimize",0,1) ;
  pc.defineInt("verbose","Verbose",0,0) ;
  pc.defineInt("doSave","Save",0,0) ;
  pc.defineInt("doTimer","Timer",0,0) ;
  pc.defineInt("plevel","PrintLevel",0,1) ;
  pc.defineInt("strat","Strategy",0,1) ;
  pc.defineInt("initHesse","InitialHesse",0,0) ;
  pc.defineInt("hesse","Hesse",0,1) ;
  pc.defineInt("minos","Minos",0,0) ;
  pc.defineInt("ext","Extended",0,2) ;
  pc.defineInt("numee","PrintEvalErrors",0,10) ;
  pc.defineInt("doWarn","Warnings",0,1) ;
  pc.defineString("mintype","Minimizer",0,"Minuit") ;
  pc.defineString("minalg","Minimizer",1,"minuit") ;
  pc.defineObject("minosSet","Minos",0,0) ;

  pc.defineMutex("FitOptions","Verbose") ;
  pc.defineMutex("FitOptions","Save") ;
  pc.defineMutex("FitOptions","Timer") ;
  pc.defineMutex("FitOptions","Strategy") ;
  pc.defineMutex("FitOptions","InitialHesse") ;
  pc.defineMutex("FitOptions","Hesse") ;
  pc.defineMutex("FitOptions","Minos") ;

  // Process and check varargs
  pc.process(cmdList) ;
  if (!pc.ok(kTRUE)) {
    return 0 ;
  }

  // Decode command line arguments
  const char* fitOpt = pc.getString("fitOpt",0,kTRUE) ;
#ifdef __ROOFIT_NOROOMINIMIZER
  const char* minType =0 ;
#else
  const char* minType = pc.getString("mintype","Minuit") ;
  const char* minAlg = pc.getString("minalg","minuit") ;
#endif
  Int_t optConst = pc.getInt("optConst") ;
  Int_t verbose  = pc.getInt("verbose") ;
  Int_t doSave   = pc.getInt("doSave") ;
  Int_t doTimer  = pc.getInt("doTimer") ;
  Int_t plevel    = pc.getInt("plevel") ;
  Int_t strat    = pc.getInt("strat") ;
  Int_t initHesse= pc.getInt("initHesse") ;
  Int_t hesse    = pc.getInt("hesse") ;
  Int_t minos    = pc.getInt("minos") ;
  Int_t numee    = pc.getInt("numee") ;
  Int_t doWarn   = pc.getInt("doWarn") ;
  const RooArgSet* minosSet = static_cast<RooArgSet*>(pc.getObject("minosSet")) ;

  RooFitResult *ret = 0 ;

#ifdef __ROOFIT_NOROOMINIMIZER
  if (true) {
#else
  if ("OldMinuit" == string(minType)) {
#endif
    // Instantiate MINUIT
    RooMinuit m(fcn) ;

    if (doWarn==0) {
      m.setNoWarn() ;
    }

    m.setPrintEvalErrors(numee) ;
    if (plevel!=1) {
      m.setPrintLevel(plevel) ;
    }

    if (optConst) {
      // Activate constant term optimization
      m.optimizeConst(optConst);
    }

    if (fitOpt) {

      // Play fit options as historically defined
      ret = m.fit(fitOpt) ;

    } else {

      if (verbose) {
   // Activate verbose options
   m.setVerbose(1) ;
      }
      if (doTimer) {
   // Activate timer options
   m.setProfile(1) ;
      }

      if (strat!=1) {
   // Modify fit strategy
   m.setStrategy(strat) ;
      }

      if (initHesse) {
   // Initialize errors with hesse
   m.hesse() ;
      }

      // Minimize using migrad
      m.migrad() ;

      if (hesse) {
   // Evaluate errors with Hesse
   m.hesse() ;
      }

      if (minos) {
   // Evaluate errs with Minos
   if (minosSet) {
     m.minos(*minosSet) ;
   } else {
     m.minos() ;
   }
      }

      // Optionally return fit result
      if (doSave) {
   string name = Form("fitresult_%s",fcn.GetName()) ;
   string title = Form("Result of fit of %s ",GetName()) ;
   ret = m.save(name.c_str(),title.c_str()) ;
      }

    }
  } else {
#ifndef __ROOFIT_NOROOMINIMIZER
    // Instantiate MINUIT
    RooMinimizer m(fcn) ;
    m.setMinimizerType(minType);

    if (doWarn==0) {
      // m.setNoWarn() ; WVE FIX THIS
    }

    m.setPrintEvalErrors(numee) ;
    if (plevel!=1) {
      m.setPrintLevel(plevel) ;
    }

    if (optConst) {
      // Activate constant term optimization
      m.optimizeConst(optConst);
    }

    if (fitOpt) {

      // Play fit options as historically defined
      ret = m.fit(fitOpt) ;

    } else {

      if (verbose) {
   // Activate verbose options
   m.setVerbose(1) ;
      }
      if (doTimer) {
   // Activate timer options
   m.setProfile(1) ;
      }

      if (strat!=1) {
   // Modify fit strategy
   m.setStrategy(strat) ;
      }

      if (initHesse) {
   // Initialize errors with hesse
   m.hesse() ;
      }

      // Minimize using migrad
      m.minimize(minType, minAlg) ;

      if (hesse) {
   // Evaluate errors with Hesse
   m.hesse() ;
      }

      if (minos) {
   // Evaluate errs with Minos
   if (minosSet) {
     m.minos(*minosSet) ;
   } else {
     m.minos() ;
   }
      }

      // Optionally return fit result
      if (doSave) {
   string name = Form("fitresult_%s",fcn.GetName()) ;
   string title = Form("Result of fit of %s ",GetName()) ;
   ret = m.save(name.c_str(),title.c_str()) ;
      }
    }
#endif
  }

  // Cleanup
  return ret ;

}


////////////////////////////////////////////////////////////////////////////////
/// Return current evaluation error logging mode.

RooAbsReal::ErrorLoggingMode RooAbsReal::evalErrorLoggingMode()
{
  return _evalErrorMode ;
}

////////////////////////////////////////////////////////////////////////////////
/// Set evaluation error logging mode. Options are
///
/// PrintErrors - Print each error through RooMsgService() as it occurs
/// CollectErrors - Accumulate errors, but do not print them. A subsequent call
///                 to printEvalErrors() will print a summary
/// CountErrors - Accumulate error count, but do not print them.
///

void RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::ErrorLoggingMode m)
{
  _evalErrorMode =  m;
}


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

void RooAbsReal::setParameterizeIntegral(const RooArgSet& paramVars)
{
  RooFIter iter = paramVars.fwdIterator() ;
  RooAbsArg* arg ;
  string plist ;
  while((arg=iter.next())) {
    if (!dependsOnValue(*arg)) {
      coutW(InputArguments) << "RooAbsReal::setParameterizeIntegral(" << GetName()
             << ") function does not depend on listed parameter " << arg->GetName() << ", ignoring" << endl ;
      continue ;
    }
    if (plist.size()>0) plist += ":" ;
    plist += arg->GetName() ;
  }
  setStringAttribute("CACHEPARAMINT",plist.c_str()) ;
}