Revision 716fb731a3c8c9b583235feddac21006994d9626 authored by Kyle Cranmer on 03 August 2011, 22:09:36 UTC, committed by Kyle Cranmer on 03 August 2011, 22:09:36 UTC
git-svn-id: http://root.cern.ch/svn/root/branches/v5-30-00-patches@40446 27541ba8-7e3a-0410-8455-c3a389f83636
1 parent 0db9ccb
RooAbsReal.cxx
/*****************************************************************************
* 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) *
*****************************************************************************/
//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// 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.
//
//
// END_HTML
//
#include <sys/types.h>
#include "RooFit.h"
#include "RooMsgService.h"
#include "RooAbsReal.h"
#include "RooAbsReal.h"
#include "RooArgSet.h"
#include "RooArgList.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"
#include "RooChi2Var.h"
#include "RooFitResult.h"
#include "RooMoment.h"
#include "RooBrentRootFinder.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 ;
RooAbsReal::ErrorLoggingMode RooAbsReal::_evalErrorMode = RooAbsReal::PrintErrors ;
Int_t RooAbsReal::_evalErrorCount = 0 ;
map<const RooAbsArg*,pair<string,list<RooAbsReal::EvalError> > > RooAbsReal::_evalErrorList ;
//_____________________________________________________________________________
RooAbsReal::RooAbsReal() : _specIntegratorConfig(0), _treeVar(kFALSE), _selectComp(kTRUE), _lastNSet(0)
{
// coverity[UNINIT_CTOR]
// Default constructor
}
//_____________________________________________________________________________
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)
{
// Constructor with unit label
setValueDirty() ;
setShapeDirty() ;
}
//_____________________________________________________________________________
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)
{
// Constructor with plot range and unit label
setValueDirty() ;
setShapeDirty() ;
}
//_____________________________________________________________________________
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), _forceNumInt(other._forceNumInt),
_treeVar(other._treeVar), _selectComp(other._selectComp), _lastNSet(0)
{
// coverity[UNINIT_CTOR]
// Copy constructor
if (other._specIntegratorConfig) {
_specIntegratorConfig = new RooNumIntConfig(*other._specIntegratorConfig) ;
} else {
_specIntegratorConfig = 0 ;
}
}
//_____________________________________________________________________________
RooAbsReal::~RooAbsReal()
{
// Destructor
if (_specIntegratorConfig) delete _specIntegratorConfig ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::operator==(Double_t value) const
{
// Equality operator comparing to a Double_t
return (getVal()==value) ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::operator==(const RooAbsArg& other)
{
// Equality operator when comparing to another RooAbsArg.
// Only functional when the other arg is a RooAbsReal
const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
return otherReal ? operator==(otherReal->getVal()) : kFALSE ;
}
//_____________________________________________________________________________
TString RooAbsReal::getTitle(Bool_t appendUnit) const
{
// Return this variable's title string. If appendUnit is true and
// this variable has units, also append a string " (<unit>)".
TString title(GetTitle());
if(appendUnit && 0 != strlen(getUnit())) {
title.Append(" (");
title.Append(getUnit());
title.Append(")");
}
return title;
}
//_____________________________________________________________________________
Double_t RooAbsReal::getVal(const RooArgSet* nset) const
{
// Return value of object. If the cache is clean, return the
// cached value, otherwise recalculate on the fly and refill
// the cache
if (nset && nset!=_lastNSet) {
((RooAbsReal*) this)->setProxyNormSet(nset) ;
_lastNSet = (RooArgSet*) nset ;
}
if (isValueDirty() || isShapeDirty()) {
_value = traceEval(nset) ;
clearValueDirty() ;
clearShapeDirty() ;
} else if (_cacheCheck) {
// Check if cache contains value that evaluate() gives now
Double_t checkValue = traceEval(nset);
if (checkValue != _value) {
// If not, print warning
coutW(Eval) << "RooAbsReal::getVal(" << GetName() << ") WARNING: cache contains " << _value
<< " but evaluate() returns " << checkValue << endl ;
// And update cache (so that we see the difference)
_value = checkValue ;
}
}
return _value ;
}
//_____________________________________________________________________________
Int_t RooAbsReal::numEvalErrorItems()
{
return _evalErrorList.size() ;
}
//_____________________________________________________________________________
RooAbsReal::EvalErrorIter RooAbsReal::evalErrorIter()
{
return _evalErrorList.begin() ;
}
//_____________________________________________________________________________
Double_t RooAbsReal::traceEval(const RooArgSet* /*nset*/) const
{
// Calculate current value of object, with error tracing wrapper
Double_t value = evaluate() ;
if (TMath::IsNaN(value)) {
logEvalError("function value is NAN") ;
}
cxcoutD(Tracing) << "RooAbsReal::getVal(" << 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 ;
}
//_____________________________________________________________________________
Int_t RooAbsReal::getAnalyticalIntegralWN(RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgSet* /*normSet*/, const char* rangeName) const
{
// 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
return _forceNumInt ? 0 : getAnalyticalIntegral(allDeps,analDeps,rangeName) ;
}
//_____________________________________________________________________________
Int_t RooAbsReal::getAnalyticalIntegral(RooArgSet& /*integSet*/, RooArgSet& /*anaIntSet*/, const char* /*rangeName*/) const
{
// 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.
return 0 ;
}
//_____________________________________________________________________________
Double_t RooAbsReal::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
{
// Implements the actual analytical integral(s) advertised by
// getAnalyticalIntegral. This functions will only be called with
// codes returned by getAnalyticalIntegral, except code zero.
// cout << "RooAbsReal::analyticalIntegralWN(" << GetName() << ") code = " << code << " normSet = " << (normSet?*normSet:RooArgSet()) << endl ;
if (code==0) return getVal(normSet) ;
return analyticalIntegral(code,rangeName) ;
}
//_____________________________________________________________________________
Double_t RooAbsReal::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
{
// Implements the actual analytical integral(s) advertised by
// getAnalyticalIntegral. This functions will only be called with
// codes returned by getAnalyticalIntegral, except code zero.
// By default no analytical integrals are implemented
coutF(Eval) << "RooAbsReal::analyticalIntegral(" << GetName() << ") code " << code << " not implemented" << endl ;
return 0 ;
}
//_____________________________________________________________________________
const char *RooAbsReal::getPlotLabel() const
{
// Get the label associated with the variable
return _label.IsNull() ? fName.Data() : _label.Data();
}
//_____________________________________________________________________________
void RooAbsReal::setPlotLabel(const char *label)
{
// Set the label associated with this variable
_label= label;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::readFromStream(istream& /*is*/, Bool_t /*compact*/, Bool_t /*verbose*/)
{
//Read object contents from stream (dummy for now)
return kFALSE ;
}
//_____________________________________________________________________________
void RooAbsReal::writeToStream(ostream& /*os*/, Bool_t /*compact*/) const
{
//Write object contents to stream (dummy for now)
}
//_____________________________________________________________________________
void RooAbsReal::printValue(ostream& os) const
{
// Print object value
os << getVal() ;
}
//_____________________________________________________________________________
void RooAbsReal::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
{
// Structure printing
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;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::isValid() const
{
// Check if current value is valid
return isValidReal(_value) ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::isValidReal(Double_t /*value*/, Bool_t /*printError*/) const
{
// Interface function to check if given value is a valid value for this object.
// This default implementation considers all values valid
return kTRUE ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createProfile(const RooArgSet& paramsOfInterest)
{
// 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
// 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) ;
}
//_____________________________________________________________________________
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
{
// 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
// 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) ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createIntegral(const RooArgSet& iset, const RooArgSet* nset,
const RooNumIntConfig* cfg, const char* rangeName) const
{
// 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
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 ;
// char* buf = new char[strlen(rangeName)+1] ;
// strlcpy(buf,rangeName) ;
// char* range = strtok(buf,",") ;
// while (range) {
// RooAbsReal* compIntegral = createIntObj(iset,nset,cfg,range) ;
// components.add(*compIntegral) ;
// range = strtok(0,",") ;
// }
// delete[] buf ;
// + ALEX
TObjArray* oa = TString(rangeName).Tokenize(",");
for( Int_t i=0; i < oa->GetEntries(); ++i) {
TObjString* os = (TObjString*) (*oa)[i];
// cout<< " ALEX:: RooAbsReal::createIntegral (" << GetName() << ") os = " << os->GetString().Data() <<endl;
if(!os) break;
RooAbsReal* compIntegral = createIntObj(iset,nset,cfg,os->GetString().Data()) ;
// cout << "just created " << compIntegral->GetName() << " for rangename = " << os->GetString().Data() << endl ;
components.add(*compIntegral) ;
}
delete oa;
// - ALEX
// cout<< " ALEX:: RooAbsReal::createIntegral (" << GetName() << ") components = " << components <<endl;
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) ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createIntObj(const RooArgSet& iset2, const RooArgSet* nset2,
const RooNumIntConfig* cfg, const char* rangeName) const
{
// Utility function for createIntegral that creates the actual integreal object
// 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 ;
}
//_____________________________________________________________________________
void RooAbsReal::findInnerMostIntegration(const RooArgSet& allObs, RooArgSet& innerObs, const char* rangeName) const
{
// 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.
// 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,kTRUE,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) ;
}
//_____________________________________________________________________________
TString RooAbsReal::integralNameSuffix(const RooArgSet& iset, const RooArgSet* nset, const char* rangeName, Bool_t omitEmpty) const
{
// Construct string with unique suffix name to give to integral object that encodes
// integrated observables, normalization observables and the integration range name
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 ;
}
//_____________________________________________________________________________
const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars,
RooArgSet*& cloneSet) const
{
// Utility function for plotOn() that creates a projection of a function or p.d.f
// to be plotted on a RooPlot.
return createPlotProjection(depVars,&projVars,cloneSet) ;
}
//_____________________________________________________________________________
const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars) const
{
// Utility function for plotOn() that creates a projection of a function or p.d.f
// to be plotted on a RooPlot.
RooArgSet* cloneSet = new RooArgSet() ;
return createPlotProjection(depVars,&projVars,cloneSet) ;
}
//_____________________________________________________________________________
const RooAbsReal *RooAbsReal::createPlotProjection(const RooArgSet &dependentVars, const RooArgSet *projectedVars,
RooArgSet *&cloneSet, const char* rangeName, const RooArgSet* condObs) const
{
// 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)
// 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= ∅
TString name = GetName() ;
name += integralNameSuffix(*projectedVars,&normSet,rangeName,kTRUE) ;
TString title(GetTitle());
title.Prepend("Projection of ");
RooRealIntegral *projected= new RooRealIntegral(name.Data(),title.Data(),*theClone,*projectedVars,&normSet,0,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;
}
// 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;
}
//_____________________________________________________________________________
TH1 *RooAbsReal::fillHistogram(TH1 *hist, const RooArgList &plotVars,
Double_t scaleFactor, const RooArgSet *projectedVars, Bool_t scaleForDensity,
const RooArgSet* condObs, Bool_t setError) const
{
// 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)
// 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),result) ;
}
//cout << "bin " << bin << " -> (" << xbin << "," << ybin << "," << zbin << ") = " << result << endl;
}
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
// cleanup
delete cloneSet;
return hist;
}
//_____________________________________________________________________________
RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSet, Double_t scaleFactor,
Bool_t correctForBinSize, Bool_t showProgress) const
{
// 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
// 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()) ;
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;
}
//_____________________________________________________________________________
TH1* RooAbsReal::createHistogram(const char* varNameList, Int_t xbins, Int_t ybins, Int_t zbins) const
{
// 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
// 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 ;
}
//_____________________________________________________________________________
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
{
// 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
//
// 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
//
// 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
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) ;
}
//_____________________________________________________________________________
TH1* RooAbsReal::createHistogram(const char *name, const RooAbsRealLValue& xvar, RooLinkedList& argList) const
{
// Internal method implementing createHistogram
// Define configuration for this method
RooCmdConfig pc(Form("RooAbsReal::createHistogram(%s)",GetName())) ;
pc.defineInt("scaling","Scaling",0,1) ;
pc.defineSet("projObs","ProjectedObservables",0,0) ;
pc.defineObject("yvar","YVar",0,0) ;
pc.defineObject("zvar","ZVar",0,0) ;
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") ;
RooLinkedList argListCreate(argList) ;
pc.stripCmdList(argListCreate,"Scaling,ProjectedObservables") ;
TH1* histo = xvar.createHistogram(name,argListCreate) ;
fillHistogram(histo,vars,1.0,intObs,doScaling,projObs,kFALSE) ;
return histo ;
}
//_____________________________________________________________________________
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
{
// 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
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) ;
}
//_____________________________________________________________________________
RooPlot* RooAbsReal::plotOn(RooPlot* frame, RooLinkedList& argList) const
{
// Internal back-end function of plotOn() with named arguments
// 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 = 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 ;
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) ;
dwa.constOptimizeTestStatistic(Activate) ;
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);
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;
}
//_____________________________________________________________________________
RooPlot* RooAbsReal::plotSliceOn(RooPlot *frame, const RooArgSet& sliceSet, Option_t* drawOptions,
Double_t scaleFactor, ScaleType stype, const RooAbsData* projData) const
{
// OBSOLETE -- RETAINED FOR BACKWARD COMPATIBILITY. Use the plotOn(frame,Slice(...)) instead
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) ;
const RooAbsReal *negProj = funcNeg->createPlotProjection(depNeg, &projectedVars, negProjCompList) ;
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;
}
//_____________________________________________________________________________
Double_t RooAbsReal::getPropagatedError(const RooFitResult& fr)
{
// 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
//
// 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) ;
}
//_____________________________________________________________________________
RooPlot* RooAbsReal::plotOnWithErrorBand(RooPlot* frame,const RooFitResult& fr, Double_t Z,const RooArgSet* params, const RooLinkedList& argList, Bool_t linMethod) const
{
// 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
RooLinkedList plotArgList(argList) ;
RooCmdConfig pc(Form("RooAbsPdf::plotOn(%s)",GetName())) ;
pc.stripCmdList(plotArgList,"VisualizeError,MoveToBack") ;
// 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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::plotSanityChecks(RooPlot* frame) const
{
// Utility function for plotOn(), perform general sanity check on frame to ensure safe plotting operations
// 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 ;
}
//_____________________________________________________________________________
void RooAbsReal::makeProjectionSet(const RooAbsArg* plotVar, const RooArgSet* allVars,
RooArgSet& projectedVars, Bool_t silent) const
{
// 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.
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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::isSelectedComp() const
{
// If true, the current pdf is a selected component (for use in plotting)
return _selectComp || _globalSelectComp ;
}
//_____________________________________________________________________________
void RooAbsReal::globalSelectComp(Bool_t flag)
{
// Global switch controlling the activation of the selectComp() functionality
_globalSelectComp = flag ;
}
//_____________________________________________________________________________
RooAbsFunc *RooAbsReal::bindVars(const RooArgSet &vars, const RooArgSet* nset, Bool_t clipInvalid) const
{
// 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 *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;
}
//_____________________________________________________________________________
void RooAbsReal::copyCache(const RooAbsArg* source, Bool_t /*valueOnly*/)
{
// 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
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("SIGNEDBYTE_TREE_BRANCH")) {
_value = other->_sbyteValue ;
} else if (source->getAttribute("UNSIGNED_INTEGER_TREE_BRANCH")) {
_value = other->_uintValue ;
}
}
setValueDirty() ;
}
//_____________________________________________________________________________
void RooAbsReal::attachToTree(TTree& t, Int_t bufSize)
{
// 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.
// 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("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 ;
}
}
//_____________________________________________________________________________
void RooAbsReal::fillTreeBranch(TTree& t)
{
// Fill the tree branch that associated with this object with its current value
// 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() ;
}
//_____________________________________________________________________________
void RooAbsReal::setTreeBranchStatus(TTree& t, Bool_t active)
{
// (De)Activate associated tree branch
TBranch* branch = t.GetBranch(cleanBranchName()) ;
if (branch) {
t.SetBranchStatus(cleanBranchName(),active?1:0) ;
}
}
//_____________________________________________________________________________
RooAbsArg *RooAbsReal::createFundamental(const char* newname) const
{
// Create a RooRealVar fundamental object with our properties. The new
// object will be created without any fit limits.
RooRealVar *fund= new RooRealVar(newname?newname:GetName(),GetTitle(),_value,getUnit());
fund->removeRange();
fund->setPlotLabel(getPlotLabel());
fund->setAttribute("fundamentalCopy");
return fund;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgProxy& a) const
{
// 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
TList nameList ;
nameList.Add(new TObjString(a.absArg()->GetName())) ;
Bool_t result = matchArgsByName(allDeps,analDeps,nameList) ;
nameList.Delete() ;
return result ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgProxy& a, const RooArgProxy& b) const
{
// 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
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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgProxy& a, const RooArgProxy& b,
const RooArgProxy& c) const
{
// 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
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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgProxy& a, const RooArgProxy& b,
const RooArgProxy& c, const RooArgProxy& d) const
{
// 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
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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
const RooArgSet& refset) const
{
// 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.
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 ;
}
//_____________________________________________________________________________
Bool_t RooAbsReal::matchArgsByName(const RooArgSet &allArgs, RooArgSet &matchedArgs,
const TList &nameList) const
{
// 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.
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;
}
//_____________________________________________________________________________
RooNumIntConfig* RooAbsReal::defaultIntegratorConfig()
{
// Returns the default numeric integration configuration for all RooAbsReals
return &RooNumIntConfig::defaultConfig() ;
}
//_____________________________________________________________________________
RooNumIntConfig* RooAbsReal::specialIntegratorConfig() const
{
// Returns the specialized integrator configuration for _this_ RooAbsReal.
// If this object has no specialized configuration, a null pointer is returned.
return _specIntegratorConfig ;
}
//_____________________________________________________________________________
RooNumIntConfig* RooAbsReal::specialIntegratorConfig(Bool_t createOnTheFly)
{
// 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
if (!_specIntegratorConfig && createOnTheFly) {
_specIntegratorConfig = new RooNumIntConfig(*defaultIntegratorConfig()) ;
}
return _specIntegratorConfig ;
}
//_____________________________________________________________________________
const RooNumIntConfig* RooAbsReal::getIntegratorConfig() const
{
// 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* config = specialIntegratorConfig() ;
if (config) return config ;
return defaultIntegratorConfig() ;
}
//_____________________________________________________________________________
RooNumIntConfig* RooAbsReal::getIntegratorConfig()
{
// 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* config = specialIntegratorConfig() ;
if (config) return config ;
return defaultIntegratorConfig() ;
}
//_____________________________________________________________________________
void RooAbsReal::setIntegratorConfig(const RooNumIntConfig& config)
{
// Set the given integrator configuration as default numeric integration
// configuration for this object
if (_specIntegratorConfig) {
delete _specIntegratorConfig ;
}
_specIntegratorConfig = new RooNumIntConfig(config) ;
}
//_____________________________________________________________________________
void RooAbsReal::setIntegratorConfig()
{
// Remove the specialized numeric integration configuration associated
// with this object
if (_specIntegratorConfig) {
delete _specIntegratorConfig ;
}
_specIntegratorConfig = 0 ;
}
//_____________________________________________________________________________
void RooAbsReal::selectNormalization(const RooArgSet*, Bool_t)
{
// 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::selectNormalizationRange(const char*, 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::setCacheCheck(Bool_t flag)
{
// Activate cache validation mode
_cacheCheck = flag ;
}
//_____________________________________________________________________________
Int_t RooAbsReal::getMaxVal(const RooArgSet& /*vars*/) const
{
// 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
return 0 ;
}
//_____________________________________________________________________________
Double_t RooAbsReal::maxVal(Int_t /*code*/) const
{
// Return maximum value for set of observables identified by code assigned
// in getMaxVal
assert(1) ;
return 0 ;
}
//_____________________________________________________________________________
void RooAbsReal::EvalError::setMessage(const char* tmp)
{
if (strlen(tmp)<1023) {
strlcpy(_msg,tmp,1023) ;
} else {
strncpy(_msg,tmp,1020);
_msg[1020]='.' ; _msg[1021]='.' ;
_msg[1022]='.' ; _msg[1023]=0 ;
}
}
//_____________________________________________________________________________
void RooAbsReal::EvalError::setServerValues(const char* tmp)
{
if (strlen(tmp)<1023) {
strlcpy(_srvval,tmp,1023) ;
} else {
strncpy(_srvval,tmp,1020);
_srvval[1020]='.' ; _srvval[1021]='.' ;
_srvval[1022]='.' ; _srvval[1023]=0 ;
}
}
//_____________________________________________________________________________
void RooAbsReal::logEvalError(const RooAbsReal* originator, const char* origName, const char* message, const char* serverValueString)
{
// Interface to insert remote error logging messages received by RooRealMPFE into current error loggin stream
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 ;
}
//_____________________________________________________________________________
void RooAbsReal::logEvalError(const char* message, const char* serverValueString) const
{
// 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.
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->name()[0]=='!') continue ;
if (first) {
first=kFALSE ;
} else {
oss << ", " ;
}
getProxy(i)->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) {
_evalErrorList[this].first = oss2.str().c_str() ;
_evalErrorList[this].second.push_back(ee) ;
}
inLogEvalError = kFALSE ;
//coutE(Tracing) << "RooAbsReal::logEvalError(" << GetName() << ") message = " << message << endl ;
}
//_____________________________________________________________________________
void RooAbsReal::clearEvalErrorLog()
{
// Clear the stack of evaluation error messages
if (_evalErrorMode==PrintErrors) {
return ;
} else if (_evalErrorMode==CollectErrors) {
_evalErrorList.clear() ;
} else {
_evalErrorCount = 0 ;
}
}
//_____________________________________________________________________________
void RooAbsReal::printEvalErrors(ostream& os, Int_t maxPerNode)
{
// 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.
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 ;
}
}
}
}
}
//_____________________________________________________________________________
Int_t RooAbsReal::numEvalErrors()
{
// Return the number of logged evaluation errors since the last clearing.
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 ;
}
//_____________________________________________________________________________
void RooAbsReal::fixAddCoefNormalization(const RooArgSet& addNormSet, Bool_t force)
{
// 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.
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 ;
}
//_____________________________________________________________________________
void RooAbsReal::fixAddCoefRange(const char* rangeName, Bool_t force)
{
// 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.
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 ;
}
//_____________________________________________________________________________
void RooAbsReal::preferredObservableScanOrder(const RooArgSet& obs, RooArgSet& orderedObs) const
{
// 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'
// Dummy implementation, do nothing
orderedObs.removeAll() ;
orderedObs.add(obs) ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createRunningIntegral(const RooArgSet& iset, const RooArgSet& nset)
{
// 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'
return createRunningIntegral(iset,RooFit::SupNormSet(nset)) ;
}
//_____________________________________________________________________________
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)
{
// 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
// 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 ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createScanRI(const RooArgSet& iset, const RooArgSet& nset, Int_t numScanBins, Int_t intOrder)
{
// Utility function for createRunningIntegral that construct an object
// implementing the numeric scanning technique for calculating the running integral
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 ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createIntRI(const RooArgSet& iset, const RooArgSet& nset)
{
// Utility function for createRunningIntegral that construct an
// object implementing the standard (analytical) integration
// technique for calculating the running integral
// 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 ;
}
//_____________________________________________________________________________
RooFunctor* RooAbsReal::functor(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
{
// Return a RooFunctor object bound to this RooAbsReal with given definition of observables
// and parameters
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) ;
}
//_____________________________________________________________________________
TF1* RooAbsReal::asTF(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
{
// Return a ROOT TF1,2,3 object bound to this RooAbsReal with given definition of observables
// and parameters
// 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(),"RooFunctor") ;
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(),"RooFunctor") ;
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(),"RooFunctor") ;
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 ;
}
//_____________________________________________________________________________
RooDerivative* RooAbsReal::derivative(RooRealVar& obs, Int_t order, Double_t eps)
{
// Return function representing first, second or third order derivative of this function
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) ;
}
//_____________________________________________________________________________
RooDerivative* RooAbsReal::derivative(RooRealVar& obs, const RooArgSet& normSet, Int_t order, Double_t eps)
{
// Return function representing first, second or third order derivative of this function
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) ;
}
//_____________________________________________________________________________
RooMoment* RooAbsReal::moment(RooRealVar& obs, Int_t order, Bool_t central, Bool_t takeRoot)
{
// Return function representing moment of function of given order. If central is
// true, the central moment is given <(x-<x>)^2>
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()) ;
return new RooMoment(name.c_str(),title.c_str(),*this,obs,order,central,takeRoot) ;
}
//_____________________________________________________________________________
RooMoment* RooAbsReal::moment(RooRealVar& obs, const RooArgSet& normObs, Int_t order, Bool_t central, Bool_t takeRoot, Bool_t intNormObs)
{
// 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.
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()) ;
return new RooMoment(name.c_str(),title.c_str(),*this,obs,normObs,order,central,takeRoot,intNormObs) ;
}
//_____________________________________________________________________________
Double_t RooAbsReal::findRoot(RooRealVar& x, Double_t xmin, Double_t xmax, Double_t yval)
{
//
// Return value of x (in range xmin,xmax) at which function equals yval.
// (Calculation is performed with Brent root finding algorithm)
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) ;
}
//_____________________________________________________________________________
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)
{
// 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.
//
//
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) ;
}
//_____________________________________________________________________________
RooFitResult* RooAbsReal::chi2FitTo(RooDataHist& data, const RooLinkedList& cmdList)
{
// Internal back-end function to steer chi2 fits
// 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 ;
}
//_____________________________________________________________________________
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)
{
// 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
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) ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createChi2(RooDataHist& data, const RooLinkedList& cmdList)
{
// Internal back-end function to create a chi2
// 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 ;
// Construct chi2
string name = Form("chi2_%s_%s",GetName(),data.GetName()) ;
return new RooChi2Var(name.c_str(),name.c_str(),*this,data,*cmds[0],*cmds[1],*cmds[2],*cmds[3],*cmds[4],*cmds[5],*cmds[6],*cmds[7]) ;
}
//_____________________________________________________________________________
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)
{
// 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.
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) ;
}
//_____________________________________________________________________________
RooFitResult* RooAbsReal::chi2FitTo(RooDataSet& xydata, const RooLinkedList& cmdList)
{
// Internal back-end function to steer chi2 fits
// 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 ;
}
//_____________________________________________________________________________
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)
{
// 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.
//
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) ;
}
//_____________________________________________________________________________
RooAbsReal* RooAbsReal::createChi2(RooDataSet& data, const RooLinkedList& cmdList)
{
// Internal back-end function to create a chi^2 from a function and a dataset
// 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) ;
}
}
//_____________________________________________________________________________
RooFitResult* RooAbsReal::chi2FitDriver(RooAbsReal& fcn, RooLinkedList& cmdList)
{
// Internal driver function for chi2 fits
// 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.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) ;
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")) ;
// 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(1) ;
}
RooFitResult *ret = 0 ;
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()) ;
}
}
// Cleanup
return ret ;
}
//_____________________________________________________________________________
RooAbsReal::ErrorLoggingMode RooAbsReal::evalErrorLoggingMode()
{
// Return current evaluation error logging mode.
return _evalErrorMode ;
}
//_____________________________________________________________________________
void RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::ErrorLoggingMode m)
{
// 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.
//
_evalErrorMode = m;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...