/***************************************************************************** * Project: RooFit * * Package: RooFitCore * * @(#)root/roofitcore:$Id: RooAbsPdf.cxx 24307 2008-06-17 13:21:58Z wouter $ * 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) * *****************************************************************************/ ////////////////////////////////////////////////////////////////////////////// // // RooAbsPdf is the abstract interface for all probability density // functions The class provides hybrid analytical/numerical // normalization for its implementations, error tracing and a MC // generator interface. // // A minimal implementation of a PDF class derived from RooAbsPdf // should overload the evaluate() function. This functions should // return PDFs value. // // // [Normalization/Integration] // // Although the normalization of a PDF is an integral part of a // probability density function, normalization is treated separately // in RooAbsPdf. The reason is that a RooAbsPdf object is more than a // PDF: it can be a building block for a more complex, composite PDF // if any of its variables are functions instead of variables. In // such cases the normalization of the composite may not be simply the // integral over the dependents of the top level PDF as these are // functions with potentially non-trivial Jacobian terms themselves. // Therefore // // --> No explicit attempt should be made to normalize // the functions output in evaluate(). // // In addition, RooAbsPdf objects do not have a static concept of what // variables are parameters and what variables are dependents (which // need to be integrated over for a correct PDF normalization). // Instead the choice of normalization is always specified each time a // normalized values is requested from the PDF via the getVal() // method. // // RooAbsPdf manages the entire normalization logic of each PDF with // help of a RooRealIntegral object, which coordinates the integration // of a given choice of normalization. By default, RooRealIntegral will // perform a fully numeric integration of all dependents. However, // PDFs can advertise one or more (partial) analytical integrals of // their function, and these will be used by RooRealIntegral, if it // determines that this is safe (i.e. no hidden Jacobian terms, // multiplication with other PDFs that have one or more dependents in // commen etc) // // To implement analytical integrals, two functions must be implemented. First, // // Int_t getAnalyticalIntegral(const RooArgSet& integSet, RooArgSet& anaIntSet) // // 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. Second, // // Double_t analyticalIntegral(Int_t code) // // Implements the actual analytical integral(s) advertised by // getAnalyticalIntegral. This functions will only be called with // codes returned by getAnalyticalIntegral, except code zero. // // The integration range for real each dependent to be integrated can // be obtained from the dependents' proxy functions min() and // max(). Never call these proxy functions for any proxy not known to // be a dependent via the integration code. Doing so may be // ill-defined, e.g. in case the proxy holds a function, and will // trigger an assert. Integrated category dependents should always be // summed over all of their states. // // // // [Direct generation of observables] // // Any PDF dependent can be generated with the accept/reject method, // but for certain PDFs more efficient methods may be implemented. To // implement direct generation of one or more observables, two // functions need to be implemented, similar to those for analytical // integrals: // // Int_t getGenerator(const RooArgSet& generateVars, RooArgSet& directVars) and // void generateEvent(Int_t code) // // The first function advertises observables that can be generated, // similar to the way analytical integrals are advertised. The second // function implements the generator for the advertised observables // // The generated dependent values should be store in the proxy // objects. For this the assignment operator can be used (i.e. xProxy // = 3.0 ). Never call assign to any proxy not known to be a dependent // via the generation code. Doing so may be ill-defined, e.g. in case // the proxy holds a function, and will trigger an assert #include "RooFit.h" #include "RooMsgService.h" #include "TClass.h" #include "Riostream.h" #include "TMath.h" #include "TObjString.h" #include "TPaveText.h" #include "TList.h" #include "TH1.h" #include "TH2.h" #include "RooAbsPdf.h" #include "RooDataSet.h" #include "RooArgSet.h" #include "RooArgProxy.h" #include "RooRealProxy.h" #include "RooRealVar.h" #include "RooGenContext.h" #include "RooPlot.h" #include "RooCurve.h" #include "RooNLLVar.h" #include "RooMinuit.h" #include "RooCategory.h" #include "RooNameReg.h" #include "RooCmdConfig.h" #include "RooGlobalFunc.h" #include "RooAddition.h" #include "RooRandom.h" #include "RooNumIntConfig.h" #include "RooProjectedPdf.h" #include "RooInt.h" #include "RooCustomizer.h" #include "RooConstraintSum.h" #include "RooParamBinning.h" #include "RooNumCdf.h" #include <string> ClassImp(RooAbsPdf) ; Int_t RooAbsPdf::_verboseEval = 0; Bool_t RooAbsPdf::_globalSelectComp = kFALSE ; Bool_t RooAbsPdf::_evalError = kFALSE ; //_____________________________________________________________________________ RooAbsPdf::RooAbsPdf() : _norm(0), _normSet(0) { // Default constructor } //_____________________________________________________________________________ RooAbsPdf::RooAbsPdf(const char *name, const char *title) : RooAbsReal(name,title), _norm(0), _normSet(0), _normMgr(this,10), _selectComp(kTRUE) { // Constructor with name and title only resetErrorCounters() ; setTraceCounter(0) ; } //_____________________________________________________________________________ RooAbsPdf::RooAbsPdf(const char *name, const char *title, Double_t plotMin, Double_t plotMax) : RooAbsReal(name,title,plotMin,plotMax), _norm(0), _normSet(0), _normMgr(this,10), _selectComp(kTRUE) { // Constructor with name, title, and plot range resetErrorCounters() ; setTraceCounter(0) ; } //_____________________________________________________________________________ RooAbsPdf::RooAbsPdf(const RooAbsPdf& other, const char* name) : RooAbsReal(other,name), _norm(0), _normSet(0), _normMgr(other._normMgr,this), _selectComp(other._selectComp) { // Copy constructor resetErrorCounters() ; setTraceCounter(other._traceCount) ; } //_____________________________________________________________________________ RooAbsPdf::~RooAbsPdf() { // Destructor } //_____________________________________________________________________________ Double_t RooAbsPdf::getVal(const RooArgSet* nset) const { // Return current value, normalizated by integrating over // the observables in 'nset'. If 'nset' is 0, the unnormalized value. // is returned. All elements of 'nset' must be lvalues // // Unnormalized values are not cached // Doing so would be complicated as _norm->getVal() could // spoil the cache and interfere with returning the cached // return value. Since unnormalized calls are typically // done in integration calls, there is no performance hit. if (!nset) { Double_t val = evaluate() ; Bool_t error = traceEvalPdf(val) ; cxcoutD(Tracing) << IsA()->GetName() << "::getVal(" << GetName() << "): value = " << val << " (unnormalized)" << endl ; if (error) { raiseEvalError() ; return 0 ; } return val ; } // Process change in last data set used Bool_t nsetChanged(kFALSE) ; if (nset!=_normSet || _norm==0) { nsetChanged = syncNormalization(nset) ; } // Return value of object. Calculated if dirty, otherwise cached value is returned. if ((isValueDirty() || nsetChanged || _norm->isValueDirty()) && operMode()!=AClean) { // cout << "doEval" << endl ; // Evaluate numerator Double_t rawVal = evaluate() ; Bool_t error = traceEvalPdf(rawVal) ; // Error checking and printing // Evaluate denominator Double_t normVal(_norm->getVal()) ; cxcoutD(Tracing) << "RooAbsPdf::getVal(" << GetName() << ") normalization integral is " << (_norm?_norm->GetName():"none") << endl ; Double_t normError(kFALSE) ; if (normVal==0.) { normError=kTRUE ; logEvalError("p.d.f normalization integral is zero") ; } // Raise global error flag if problems occur if (normError||error) raiseEvalError() ; _value = normError ? 0 : (rawVal / normVal) ; cxcoutD(Tracing) << "RooAbsPdf::getVal(" << GetName() << ") recalculating, new value = " << rawVal << "/" << normVal << " = " << _value << endl ; clearValueDirty() ; //setValueDirty(kFALSE) ; clearShapeDirty() ; //setShapeDirty(kFALSE) ; } if (_traceCount>0) { cxcoutD(Tracing) << "[" << _traceCount << "] " ; Int_t tmp = _traceCount ; _traceCount = 0 ; printStream(ccoutD(Tracing),kName|kValue|kArgs,kSingleLine) ; _traceCount = tmp-1 ; } return _value ; } //_____________________________________________________________________________ Double_t RooAbsPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const { // Analytical integral with normalization (see RooAbsReal::analyticalIntegralWN() for further information) // // This function applies the normalization specified by 'normSet' to the integral returned // by RooAbsReal::analyticalIntegral(). The passthrough scenario (code=0) is also changed // to return a normalized answer cxcoutD(Eval) << "RooAbsPdf::analyticalIntegralWN(" << GetName() << ") code = " << code << " normset = " << (normSet?*normSet:RooArgSet()) << endl ; if (code==0) return getVal(normSet) ; if (normSet) { return analyticalIntegral(code,rangeName) / getNorm(normSet) ; } else { return analyticalIntegral(code,rangeName) ; } } //_____________________________________________________________________________ Bool_t RooAbsPdf::traceEvalPdf(Double_t value) const { // Check that passed value is positive and not 'not-a-number'. If // not, print an error, until the error counter reaches its set // maximum. cxcoutD(Tracing) << "RooAbsPdf::traceEvalPdf(" << GetName() << ") value = " << value << endl ; // check for a math error or negative value Bool_t error= isnan(value) || (value < 0); if (isnan(value)) { logEvalError(Form("p.d.f value is Not-a-Number (%f), forcing value to zero",value)) ; } if (value<0) { logEvalError(Form("p.d.f value is less than zero (%f), forcing value to zero",value)) ; } // do nothing if we are no longer tracing evaluations and there was no error if(!error) return error ; // otherwise, print out this evaluations input values and result if(++_errorCount <= 10) { cxcoutD(Tracing) << "*** Evaluation Error " << _errorCount << " "; if(_errorCount == 10) cxcoutD(Tracing) << "(no more will be printed) "; } else { return error ; } Print() ; return error ; } //_____________________________________________________________________________ Double_t RooAbsPdf::getNorm(const RooArgSet* nset) const { // Return the integral of this PDF over all observables listed in 'nset'. if (!nset) return 1 ; syncNormalization(nset,kTRUE) ; if (_verboseEval>1) cxcoutD(Tracing) << IsA()->GetName() << "::getNorm(" << GetName() << "): norm(" << _norm << ") = " << _norm->getVal() << endl ; Double_t ret = _norm->getVal() ; if (ret==0.) { if(++_errorCount <= 10) { coutW(Eval) << "RooAbsPdf::getNorm(" << GetName() << ":: WARNING normalization is zero, nset = " ; nset->Print("1") ; if(_errorCount == 10) coutW(Eval) << "RooAbsPdf::getNorm(" << GetName() << ") INFO: no more messages will be printed " << endl ; } } return ret ; } //_____________________________________________________________________________ const RooAbsReal* RooAbsPdf::getNormObj(const RooArgSet* nset, const RooArgSet* iset, const TNamed* rangeName) const { // Return pointer to RooAbsReal object that implements calculation of integral over observables iset in range // rangeName, optionally taking the integrand normalized over observables nset // Check normalization is already stored CacheElem* cache = (CacheElem*) _normMgr.getObj(nset,iset,0,rangeName) ; if (cache) { return cache->_norm ; } // If not create it now RooArgSet* depList = getObservables(iset) ; RooAbsReal* norm = createIntegral(*depList,*nset, *getIntegratorConfig(), RooNameReg::str(rangeName)) ; delete depList ; // Store it in the cache cache = new CacheElem(*norm) ; _normMgr.setObj(nset,iset,cache,rangeName) ; // And return the newly created integral return norm ; } //_____________________________________________________________________________ Bool_t RooAbsPdf::syncNormalization(const RooArgSet* nset, Bool_t adjustProxies) const { // Verify that the normalization integral cached with this PDF // is valid for given set of normalization observables // // If not, the cached normalization integral (if any) is deleted // and a new integral is constructed for use with 'nset' // Elements in 'nset' can be discrete and real, but must be lvalues // // For functions that declare to be self-normalized by overloading the // selfNormalized() function, a unit normalization is always constructed _normSet = (RooArgSet*) nset ; // Check if data sets are identical CacheElem* cache = (CacheElem*) _normMgr.getObj(nset) ; if (cache) { Bool_t nsetChanged = (_norm!=cache->_norm) ; _norm = cache->_norm ; if (nsetChanged && adjustProxies) { // Update dataset pointers of proxies ((RooAbsPdf*) this)->setProxyNormSet(nset) ; } return nsetChanged ; } // Update dataset pointers of proxies if (adjustProxies) { ((RooAbsPdf*) this)->setProxyNormSet(nset) ; } RooArgSet* depList = getObservables(nset) ; if (_verboseEval>0) { if (!selfNormalized()) { cxcoutD(Tracing) << IsA()->GetName() << "::syncNormalization(" << GetName() << ") recreating normalization integral " << endl ; if (depList) depList->printStream(ccoutD(Tracing),kName|kValue|kArgs,kSingleLine) ; else ccoutD(Tracing) << "<none>" << endl ; } else { cxcoutD(Tracing) << IsA()->GetName() << "::syncNormalization(" << GetName() << ") selfNormalized, creating unit norm" << endl; } } // Destroy old normalization & create new if (selfNormalized() || !dependsOn(*depList)) { TString ntitle(GetTitle()) ; ntitle.Append(" Unit Normalization") ; TString nname(GetName()) ; nname.Append("_UnitNorm") ; _norm = new RooRealVar(nname.Data(),ntitle.Data(),1) ; } else { _norm = createIntegral(*depList,*getIntegratorConfig()) ; } // Register new normalization with manager (takes ownership) cache = new CacheElem(*_norm) ; _normMgr.setObj(nset,cache) ; delete depList ; return kTRUE ; } //_____________________________________________________________________________ Bool_t RooAbsPdf::traceEvalHook(Double_t value) const { // WVE 08/21/01 Probably obsolete now. // Floating point error checking and tracing for given float value // check for a math error or negative value Bool_t error= isnan(value) || (value < 0); // do nothing if we are no longer tracing evaluations and there was no error if(!error && _traceCount <= 0) return error ; // otherwise, print out this evaluations input values and result if(error && ++_errorCount <= 10) { cxcoutD(Tracing) << "*** Evaluation Error " << _errorCount << " "; if(_errorCount == 10) ccoutD(Tracing) << "(no more will be printed) "; } else if(_traceCount > 0) { ccoutD(Tracing) << '[' << _traceCount-- << "] "; } else { return error ; } Print() ; return error ; } //_____________________________________________________________________________ void RooAbsPdf::resetErrorCounters(Int_t resetValue) { // Reset error counter to given value, limiting the number // of future error messages for this pdf to 'resetValue' _errorCount = resetValue ; _negCount = resetValue ; } //_____________________________________________________________________________ void RooAbsPdf::setTraceCounter(Int_t value, Bool_t allNodes) { // Reset trace counter to given value, limiting the // number of future trace messages for this pdf to 'value' if (!allNodes) { _traceCount = value ; return ; } else { RooArgList branchList ; branchNodeServerList(&branchList) ; TIterator* iter = branchList.createIterator() ; RooAbsArg* arg ; while((arg=(RooAbsArg*)iter->Next())) { RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ; if (pdf) pdf->setTraceCounter(value,kFALSE) ; } delete iter ; } } //_____________________________________________________________________________ Double_t RooAbsPdf::getLogVal(const RooArgSet* nset) const { // Return the log of the current value with given normalization // An error message is printed if the argument of the log is negative. Double_t prob = getVal(nset) ; if(prob <= 0) { logEvalError("getLogVal() top-level p.d.f evaluates to zero or negative number") ; return 0; } return log(prob); } //_____________________________________________________________________________ Double_t RooAbsPdf::extendedTerm(UInt_t observed, const RooArgSet* nset) const { // Returned the extended likelihood term (Nexpect - Nobserved*log(NExpected) // of this PDF for the given number of observed events // // For successfull operation the PDF implementation must indicate // it is extendable by overloading canBeExtended() and must // implemented the expectedEvents() function. // check if this PDF supports extended maximum likelihood fits if(!canBeExtended()) { coutE(InputArguments) << fName << ": this PDF does not support extended maximum likelihood" << endl; return 0; } Double_t expected= expectedEvents(nset); if(expected < 0) { coutE(InputArguments) << fName << ": calculated negative expected events: " << expected << endl; return 0; } // calculate and return the negative log-likelihood of the Poisson // factor for this dataset, dropping the constant log(observed!) Double_t extra= expected - observed*log(expected); Bool_t trace(kFALSE) ; if(trace) { cxcoutD(Tracing) << fName << "::extendedTerm: expected " << expected << " events, got " << observed << " events. extendedTerm = " << extra << endl; } return extra; } //_____________________________________________________________________________ RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, RooCmdArg arg1, RooCmdArg arg2, RooCmdArg arg3, RooCmdArg arg4, RooCmdArg arg5, RooCmdArg arg6, RooCmdArg arg7, RooCmdArg arg8) { // Fit PDF to given dataset. If dataset is unbinned, an unbinned maximum likelihood is performed. If the dataset // is binned, a binned maximum likelihood is performed. By default the fit is executed through the MINUIT // commands MIGRAD, HESSE and MINOS in succession. // // The following named arguments are supported // // Options to control construction of -log(L) // ------------------------------------------ // ConditionalObservables(const RooArgSet& set) -- Do not normalize PDF over listed observables // Extended(Bool_t flag) -- Add extended likelihood term, off by default // 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. // SumCoefRange(const char* name) -- Set the range in which to interpret the coefficients of RooAddPdf components // NumCPU(int num) -- Parallelize NLL calculation on num CPUs // Optimize(Bool_t flag) -- Activate constant term optimization (on by default) // SplitRange(Bool_t flag) -- Use separate fit ranges in a simultaneous fit. Actual range name for each // subsample is assumed to by rangeName_{indexState} where indexState // is the state of the master index category of the simultaneous fit // Contrain(const RooArgSet&pars) -- Include constraints to listed parameters in likelihood using internal constrains in p.d.f // ExternalConstraints(const RooArgSet& ) -- Include given external constraints to likelihood // // 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 // 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 fitTo(data,l) ; } //_____________________________________________________________________________ RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, const RooLinkedList& cmdList) { // Fit PDF to given dataset. If dataset is unbinned, an unbinned maximum likelihood is performed. If the dataset // is binned, a binned maximum likelihood is performed. By default the fit is executed through the MINUIT // commands MIGRAD, HESSE and MINOS in succession. // // See RooAbsPdf::fitTo(RooAbsData& data, RooCmdArg arg1, RooCmdArg arg2, RooCmdArg arg3, RooCmdArg arg4, // RooCmdArg arg5, RooCmdArg arg6, RooCmdArg arg7, RooCmdArg arg8) // // for documentation of options // Select the pdf-specific commands RooCmdConfig pc(Form("RooAbsPdf::fitTo(%s)",GetName())) ; pc.defineString("fitOpt","FitOptions",0,"") ; pc.defineString("rangeName","RangeWithName",0,"",kTRUE) ; pc.defineString("addCoefRange","SumCoefRange",0,"") ; pc.defineDouble("rangeLo","Range",0,-999.) ; pc.defineDouble("rangeHi","Range",1,-999.) ; pc.defineInt("splitRange","SplitRange",0,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,1) ; pc.defineInt("ext","Extended",0,0) ; pc.defineInt("numcpu","NumCPU",0,1) ; pc.defineInt("numee","PrintEvalErrors",0,10) ; pc.defineObject("projDepSet","ProjectedObservables",0,0) ; pc.defineObject("minosSet","Minos",0,0) ; pc.defineObject("cPars","Constrain",0,0) ; pc.defineObject("extCons","ExternalConstraints",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") ; pc.defineMutex("Range","RangeWithName") ; // 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) ; const char* rangeName = pc.getString("rangeName",0,kTRUE) ; const char* addCoefRangeName = pc.getString("addCoefRange",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 ext = pc.getInt("ext") ; Int_t numcpu = pc.getInt("numcpu") ; Int_t splitr = pc.getInt("splitRange") ; Int_t numee = pc.getInt("numee") ; const RooArgSet* minosSet = static_cast<RooArgSet*>(pc.getObject("minosSet")) ; const RooArgSet* cPars = static_cast<RooArgSet*>(pc.getObject("cPars")) ; const RooArgSet* extCons = static_cast<RooArgSet*>(pc.getObject("extCons")) ; if (pc.hasProcessed("Range")) { Double_t rangeLo = pc.getDouble("rangeLo") ; Double_t rangeHi = pc.getDouble("rangeHi") ; // Create range with name 'fit' with above limits on all observables RooArgSet* obs = getObservables(&data) ; TIterator* iter = obs->createIterator() ; RooAbsArg* arg ; while((arg=(RooAbsArg*)iter->Next())) { RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ; if (rrv) rrv->setRange("fit",rangeLo,rangeHi) ; } // Set range name to be fitted to "fit" rangeName = "fit" ; } RooArgSet projDeps ; RooArgSet* tmp = (RooArgSet*) pc.getObject("projDepSet") ; if (tmp) { projDeps.add(*tmp) ; } // Construct NLL RooAbsReal* nll ; if (!rangeName || strchr(rangeName,',')==0) { // Simple case: default range, or single restricted range nll = new RooNLLVar("nll","-log(likelihood)",*this,data,projDeps,ext,rangeName,addCoefRangeName,numcpu,kFALSE,plevel!=-1,splitr) ; } else { // Composite case: multiple ranges RooArgList nllList ; char* buf = new char[strlen(rangeName)+1] ; strcpy(buf,rangeName) ; char* token = strtok(buf,",") ; while(token) { RooAbsReal* nllComp = new RooNLLVar(Form("nll_%s",token),"-log(likelihood)",*this,data,projDeps,ext,token,addCoefRangeName,numcpu,kFALSE,plevel!=-1,splitr) ; nllList.add(*nllComp) ; token = strtok(0,",") ; } delete[] buf ; nll = new RooAddition("nll","-log(likelihood)",nllList,kTRUE) ; } // Collect internal and external constraint specifications RooArgSet allConstraints ; if (cPars) { RooArgSet* constraints = getAllConstraints(*data.get(),*cPars) ; allConstraints.add(*constraints) ; delete constraints ; } if (extCons) { allConstraints.add(*extCons) ; } // Include constraints, if any, in likelihood RooAbsReal* nllCons ; if (allConstraints.getSize()>0) { coutI(Minimization) << " Including the following contraint terms in minimization: " << allConstraints << endl ; nllCons = new RooConstraintSum("nllCons","nllCons",allConstraints) ; RooAbsReal* orignll = nll ; nll = new RooAddition("nllWithCons","nllWithCons",RooArgSet(*nll,*nllCons)) ; nll->addOwnedComponents(*orignll) ; } // Instantiate MINUIT RooMinuit m(*nll) ; 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) { ret = m.save() ; } } // Cleanup delete nll ; return ret ; } //_____________________________________________________________________________ RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, Option_t *fitOpt, Option_t *optOpt, const char* fitRange) { // OLD STYLE INTERFACE, PLEASE USE NEW INTERFACE fitTo(RooAbsData& data, RooCmdArg arg1,...,RooCmdArg arg8) return fitTo(data,RooArgSet(),fitOpt,optOpt,fitRange) ; } //_____________________________________________________________________________ RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, const RooArgSet& projDeps, Option_t *fitOpt, Option_t *optOpt, const char* fitRange) { // Fit this PDF to given data set // // OLD STYLE INTERFACE, PLEASE USE NEW INTERFACE fitTo(RooAbsData& data, RooCmdArg arg1,...,RooCmdArg arg8) // // The dataset can be either binned, in which case a binned maximum likelihood fit // is performed, or unbinned, in which case an unbinned maximum likelihood fit is performed // // Available fit options: // // "m" = MIGRAD only, i.e. no MINOS // "s" = estimate step size with HESSE before starting MIGRAD // "h" = run HESSE after MIGRAD // "e" = Perform extended MLL fit // "0" = Run MIGRAD with strategy MINUIT 0 (no correlation matrix calculation at end) // Does not apply to HESSE or MINOS, if run afterwards. // // "q" = Switch off verbose mode // "l" = Save log file with parameter values at each MINUIT step // "v" = Show changed parameters at each MINUIT step // "t" = Time fit // "r" = Save fit output in RooFitResult object (return value is object RFR pointer) // // Available optimizer options // // "c" = Cache and precalculate components of PDF that exclusively depend on constant parameters // "2" = Do NLL calculation in multi-processor mode on 2 processors // "3" = Do NLL calculation in multi-processor mode on 3 processors // "4" = Do NLL calculation in multi-processor mode on 4 processors // // The actual fit is performed to a temporary copy of both PDF and data set. Several optimization // algorithm are run to increase the efficiency of the likelihood calculation and may increase // the speed of complex fits up to an order of magnitude. All optimizations are exact, i.e the fit result // of any fit should _exactly_ the same with and without optimization. We strongly encourage // to stick to the default optimizer setting (all on). If for any reason you see a difference in the result // with and without optimizer, please file a bug report. // // The function always return null unless the "r" fit option is specified. In that case a pointer to a RooFitResult // is returned. The RooFitResult object contains the full fit output, including the correlation matrix. // Parse option strings TString fopt(fitOpt) ; TString oopt(optOpt) ; fopt.ToLower() ; oopt.ToLower() ; Bool_t extended = fopt.Contains("e") ; // Bool_t saveRes = fopt.Contains("r") ; Bool_t cOpt = oopt.Contains("p") || // for backward compatibility oopt.Contains("c") ; Bool_t blindfit = fopt.Contains("b") ; Int_t ncpu = 1 ; if (oopt.Contains("2")) ncpu=2 ; if (oopt.Contains("3")) ncpu=3 ; if (oopt.Contains("4")) ncpu=4 ; if (oopt.Contains("5")) ncpu=5 ; if (oopt.Contains("6")) ncpu=6 ; if (oopt.Contains("7")) ncpu=7 ; if (oopt.Contains("8")) ncpu=8 ; if (oopt.Contains("9")) ncpu=9 ; // Construct NLL RooNLLVar nll("nll","-log(likelihood)",*this,data,projDeps,extended,fitRange,0,ncpu) ; // Minimize NLL RooMinuit m(nll) ; if(blindfit) m.setPrintLevel(-1); if (cOpt) m.optimizeConst(1) ; return m.fit(fopt) ; } //_____________________________________________________________________________ void RooAbsPdf::printValue(ostream& os) const { // Print value of p.d.f, also print normalization integral that was last used, if any if (_norm) { os << evaluate() << "/" << _norm->getVal() ; } else { os << evaluate() ; } } //_____________________________________________________________________________ void RooAbsPdf::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const { // Print multi line detailed information of this RooAbsPdf RooAbsReal::printMultiline(os,contents,verbose,indent); os << indent << "--- RooAbsPdf ---" << endl; os << indent << "Cached value = " << _value << endl ; if (_norm) { os << indent << " Normalization integral: " << endl ; TString moreIndent(indent) ; moreIndent.Append(" ") ; _norm->printStream(os,kName|kAddress|kTitle|kValue|kArgs,kSingleLine,moreIndent.Data()) ; } } //_____________________________________________________________________________ RooAbsGenContext* RooAbsPdf::genContext(const RooArgSet &vars, const RooDataSet *prototype, const RooArgSet* auxProto, Bool_t verbose) const { // Interface function to create a generator context from a p.d.f. This default // implementation returns a 'standard' context that works for any p.d.f return new RooGenContext(*this,vars,prototype,auxProto,verbose) ; } //_____________________________________________________________________________ RooDataSet *RooAbsPdf::generate(const RooArgSet& whatVars, Int_t nEvents, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3,const RooCmdArg& arg4, const RooCmdArg& arg5) { // Generate a new dataset containing the specified variables with events sampled from our distribution. // Generate the specified number of events or expectedEvents() if not specified. // // Any variables of this PDF that are not in whatVars will use their // current values and be treated as fixed parameters. Returns zero // in case of an error. The caller takes ownership of the returned // dataset. // // The following named arguments are supported // // Verbose(Bool_t flag) -- Print informational messages during event generation // Extended() -- The actual number of events generated will be sampled from a Poisson distribution // with mu=nevt. For use with extended maximum likelihood fits // ProtoData(const RooDataSet& data, -- Use specified dataset as prototype dataset. If randOrder is set to true // Bool_t randOrder) the order of the events in the dataset will be read in a random order // if the requested number of events to be generated does not match the // number of events in the prototype dataset // // If ProtoData() is used, the specified existing dataset as a prototype: the new dataset will contain // the same number of events as the prototype (unless otherwise specified), and any prototype variables not in // whatVars will be copied into the new dataset for each generated event and also used to set our PDF parameters. // The user can specify a number of events to generate that will override the default. The result is a // copy of the prototype dataset with only variables in whatVars randomized. Variables in whatVars that // are not in the prototype will be added as new columns to the generated dataset. return generate(whatVars,RooFit::NumEvents(nEvents),arg1,arg2,arg3,arg4,arg5) ; } //_____________________________________________________________________________ RooDataSet *RooAbsPdf::generate(const RooArgSet& whatVars, const RooCmdArg& arg1,const RooCmdArg& arg2, const RooCmdArg& arg3,const RooCmdArg& arg4, const RooCmdArg& arg5,const RooCmdArg& arg6) { // Generate a new dataset containing the specified variables with events sampled from our distribution. // Generate the specified number of events or expectedEvents() if not specified. // // Any variables of this PDF that are not in whatVars will use their // current values and be treated as fixed parameters. Returns zero // in case of an error. The caller takes ownership of the returned // dataset. // // The following named arguments are supported // // Name(const char* name) -- Name of the output dataset // Verbose(Bool_t flag) -- Print informational messages during event generation // NumEvent(int nevt) -- Generate specified number of events // Extended() -- The actual number of events generated will be sampled from a Poisson distribution // with mu=nevt. For use with extended maximum likelihood fits // ProtoData(const RooDataSet& data, -- Use specified dataset as prototype dataset. If randOrder is set to true // Bool_t randOrder, the order of the events in the dataset will be read in a random order // Bool_t resample) if the requested number of events to be generated does not match the // number of events in the prototype dataset. If resample is also set to // true, the prototype dataset will be resampled rather than be strictly // reshuffled. In this mode events of the protodata may be used more than // once. // // If ProtoData() is used, the specified existing dataset as a prototype: the new dataset will contain // the same number of events as the prototype (unless otherwise specified), and any prototype variables not in // whatVars will be copied into the new dataset for each generated event and also used to set our PDF parameters. // The user can specify a number of events to generate that will override the default. The result is a // copy of the prototype dataset with only variables in whatVars randomized. Variables in whatVars that // are not in the prototype will be added as new columns to the generated dataset. // Select the pdf-specific commands RooCmdConfig pc(Form("RooAbsPdf::generate(%s)",GetName())) ; pc.defineObject("proto","PrototypeData",0,0) ; pc.defineString("dsetName","Name",0,"") ; pc.defineInt("randProto","PrototypeData",0,0) ; pc.defineInt("resampleProto","PrototypeData",1,0) ; pc.defineInt("verbose","Verbose",0,0) ; pc.defineInt("extended","Extended",0,0) ; pc.defineInt("nEvents","NumEvents",0,0) ; // Process and check varargs pc.process(arg1,arg2,arg3,arg4,arg5,arg6) ; if (!pc.ok(kTRUE)) { return 0 ; } // Decode command line arguments RooDataSet* protoData = static_cast<RooDataSet*>(pc.getObject("proto",0)) ; const char* dsetName = pc.getString("dsetName") ; Int_t nEvents = pc.getInt("nEvents") ; Bool_t verbose = pc.getInt("verbose") ; Bool_t randProto = pc.getInt("randProto") ; Bool_t resampleProto = pc.getInt("resampleProto") ; Bool_t extended = pc.getInt("extended") ; if (extended) { nEvents = RooRandom::randomGenerator()->Poisson(nEvents==0?expectedEvents(&whatVars):nEvents) ; cxcoutI(Generation) << " Extended mode active, number of events generated (" << nEvents << ") is Poisson fluctuation on " << GetName() << "::expectedEvents() = " << nEvents << endl ; } else if (nEvents==0) { cxcoutI(Generation) << "No number of events specified , number of events generated is " << GetName() << "::expectedEvents() = " << expectedEvents(&whatVars)<< endl ; } if (extended && protoData && !randProto) { cxcoutI(Generation) << "WARNING Using generator option Extended() (Poisson distribution of #events) together " << "with a prototype dataset implies incomplete sampling or oversampling of proto data. " << "Set randomize flag in ProtoData() option to randomize prototype dataset order and thus " << "to randomize the set of over/undersampled prototype events for each generation cycle." << endl ; } // Forward to appropiate implementation RooDataSet* data ; if (protoData) { data = generate(whatVars,*protoData,nEvents,verbose,randProto,resampleProto) ; } else { data = generate(whatVars,nEvents,verbose) ; } // Rename dataset to given name if supplied if (dsetName && strlen(dsetName)>0) { data->SetName(dsetName) ; } return data ; } //_____________________________________________________________________________ RooDataSet *RooAbsPdf::generate(const RooArgSet &whatVars, Int_t nEvents, Bool_t verbose) const { // Generate a new dataset containing the specified variables with // events sampled from our distribution. Generate the specified // number of events or else try to use expectedEvents() if nEvents <= 0. // Any variables of this PDF that are not in whatVars will use their // current values and be treated as fixed parameters. Returns zero // in case of an error. The caller takes ownership of the returned // dataset. RooDataSet *generated = 0; RooAbsGenContext *context= genContext(whatVars,0,0,verbose); if(0 != context && context->isValid()) { generated= context->generate(nEvents); } else { coutE(Generation) << "RooAbsPdf::generate(" << GetName() << ") cannot create a valid context" << endl; } if(0 != context) delete context; return generated; } //_____________________________________________________________________________ RooDataSet *RooAbsPdf::generate(const RooArgSet &whatVars, const RooDataSet &prototype, Int_t nEvents, Bool_t verbose, Bool_t randProtoOrder, Bool_t resampleProto) const { // Generate a new dataset with values of the whatVars variables // sampled from our distribution. Use the specified existing dataset // as a prototype: the new dataset will contain the same number of // events as the prototype (by default), and any prototype variables not in // whatVars will be copied into the new dataset for each generated // event and also used to set our PDF parameters. The user can specify a // number of events to generate that will override the default. The result is a // copy of the prototype dataset with only variables in whatVars // randomized. Variables in whatVars that are not in the prototype // will be added as new columns to the generated dataset. Returns // zero in case of an error. The caller takes ownership of the // returned dataset. RooDataSet *generated = 0; RooAbsGenContext *context= genContext(whatVars,&prototype,0,verbose); // Resampling implies reshuffling in the implementation if (resampleProto) { randProtoOrder=kTRUE ; } if (randProtoOrder && prototype.numEntries()!=nEvents) { coutI(Generation) << "RooAbsPdf::generate (Re)randomizing event order in prototype dataset (Nevt=" << nEvents << ")" << endl ; Int_t* newOrder = randomizeProtoOrder(prototype.numEntries(),nEvents,resampleProto) ; context->setProtoDataOrder(newOrder) ; delete[] newOrder ; } if(0 != context && context->isValid()) { generated= context->generate(nEvents); } else { coutE(Generation) << "RooAbsPdf::generate(" << GetName() << ") cannot create a valid context" << endl; } if(0 != context) delete context; return generated; } //_____________________________________________________________________________ Int_t* RooAbsPdf::randomizeProtoOrder(Int_t nProto, Int_t, Bool_t resampleProto) const { // Return lookup table with randomized access order for prototype events, // given nProto prototype data events and nGen events that will actually // be accessed // Make unsorted linked list of indeces RooLinkedList l ; Int_t i ; for (i=0 ; i<nProto ; i++) { l.Add(new RooInt(i)) ; } // Make output list Int_t* lut = new Int_t[nProto] ; // Randomly samply input list into output list if (!resampleProto) { // In this mode, randomization is a strict reshuffle of the order for (i=0 ; i<nProto ; i++) { Int_t iran = RooRandom::integer(nProto-i) ; RooInt* sample = (RooInt*) l.At(iran) ; lut[i] = *sample ; l.Remove(sample) ; delete sample ; } } else { // In this mode, we resample, i.e. events can be used more than once for (i=0 ; i<nProto ; i++) { lut[i] = RooRandom::integer(nProto); } } return lut ; } //_____________________________________________________________________________ Int_t RooAbsPdf::getGenerator(const RooArgSet &/*directVars*/, RooArgSet &/*generatedVars*/, Bool_t /*staticInitOK*/) const { // Load generatedVars with the subset of directVars that we can generate events for, // and return a code that specifies the generator algorithm we will use. A code of // zero indicates that we cannot generate any of the directVars (in this case, nothing // should be added to generatedVars). Any non-zero codes will be passed to our generateEvent() // implementation, but otherwise its value is arbitrary. The default implemetation of // this method returns zero. Subclasses will usually implement this method using the // matchArgs() methods to advertise the algorithms they provide. return 0 ; } //_____________________________________________________________________________ void RooAbsPdf::initGenerator(Int_t /*code*/) { // Interface for one-time initialization to setup the generator for the specified code. } //_____________________________________________________________________________ void RooAbsPdf::generateEvent(Int_t /*code*/) { // Interface for generation of anan event using the algorithm // corresponding to the specified code. The meaning of each code is // defined by the getGenerator() implementation. The default // implementation does nothing. } //_____________________________________________________________________________ Bool_t RooAbsPdf::isDirectGenSafe(const RooAbsArg& arg) const { // Check if given observable can be safely generated using the // pdfs internal generator mechanism (if that existsP). Observables // on which a PDF depends via more than route are not safe // for use with internal generators because they introduce // correlations not known to the internal generator // Arg must be direct server of self if (!findServer(arg.GetName())) return kFALSE ; // There must be no other dependency routes TIterator* sIter = serverIterator() ; const RooAbsArg *server = 0; while((server=(const RooAbsArg*)sIter->Next())) { if(server == &arg) continue; if(server->dependsOn(arg)) { delete sIter ; return kFALSE ; } } delete sIter ; return kTRUE ; } //_____________________________________________________________________________ RooPlot* RooAbsPdf::plotOn(RooPlot* frame, RooLinkedList& cmdList) 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 // 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 // // Misc content control // -------------------- // 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 // ---------------- // 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 // Invisble(Bool_t flag) -- Add curve to frame, but do not display. Useful in combination AddTo() // Sanity checks if (plotSanityChecks(frame)) return frame ; // Select the pdf-specific commands RooCmdConfig pc(Form("RooAbsPdf::plotOn(%s)",GetName())) ; pc.defineDouble("scaleFactor","Normalization",0,1.0) ; pc.defineInt("scaleType","Normalization",0,RooAbsPdf::Relative) ; pc.defineObject("compSet","SelectCompSet",0) ; pc.defineString("compSpec","SelectCompSpec",0) ; pc.defineObject("asymCat","Asymmetry",0) ; pc.defineDouble("rangeLo","Range",0,-999.) ; pc.defineDouble("rangeHi","Range",1,-999.) ; pc.defineString("rangeName","RangeWithName",0,"") ; pc.defineInt("rangeAdjustNorm","Range",0,0) ; pc.defineInt("rangeWNAdjustNorm","RangeWithName",0,0) ; pc.defineMutex("SelectCompSet","SelectCompSpec") ; pc.defineMutex("Range","RangeWithName") ; pc.allowUndefined() ; // unknowns may be handled by RooAbsReal // Process and check varargs pc.process(cmdList) ; if (!pc.ok(kTRUE)) { return frame ; } // Decode command line arguments ScaleType stype = (ScaleType) pc.getInt("scaleType") ; Double_t scaleFactor = pc.getDouble("scaleFactor") ; const RooAbsCategoryLValue* asymCat = (const RooAbsCategoryLValue*) pc.getObject("asymCat") ; const char* compSpec = pc.getString("compSpec") ; const RooArgSet* compSet = (const RooArgSet*) pc.getObject("compSet") ; Bool_t haveCompSel = (strlen(compSpec)>0 || compSet) ; // Suffix for curve name TString nameSuffix ; if (compSpec && strlen(compSpec)>0) { nameSuffix.Append("_Comp[") ; nameSuffix.Append(compSpec) ; nameSuffix.Append("]") ; } else if (compSet) { nameSuffix.Append("_Comp[") ; nameSuffix.Append(compSet->contentsString().c_str()) ; nameSuffix.Append("]") ; } // Remove PDF-only commands from command list pc.stripCmdList(cmdList,"SelectCompSet,SelectCompSpec") ; // Adjust normalization, if so requested if (asymCat) { RooCmdArg cnsuffix("CurveNameSuffix",0,0,0,0,nameSuffix.Data(),0,0,0) ; cmdList.Add(&cnsuffix); return RooAbsReal::plotOn(frame,cmdList) ; } // More sanity checks Double_t nExpected(1) ; if (stype==RelativeExpected) { if (!canBeExtended()) { coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName() << "): ERROR the 'Expected' scale option can only be used on extendable PDFs" << endl ; return frame ; } nExpected = expectedEvents(frame->getNormVars()) ; } if (stype != Raw) { if (frame->getFitRangeNEvt() && stype==Relative) { Bool_t hasCustomRange(kFALSE), adjustNorm(kFALSE) ; Double_t rangeLo(0), rangeHi(0) ; // Retrieve plot range to be able to adjust normalization to data if (pc.hasProcessed("Range")) { rangeLo = pc.getDouble("rangeLo") ; rangeHi = pc.getDouble("rangeHi") ; adjustNorm = pc.getInt("rangeAdjustNorm") ; hasCustomRange = kTRUE ; } else if (pc.hasProcessed("RangeWithName")) { rangeLo = frame->getPlotVar()->getMin(pc.getString("rangeName",0,kTRUE)) ; rangeHi = frame->getPlotVar()->getMax(pc.getString("rangeName",0,kTRUE)) ; adjustNorm = pc.getInt("rangeWNAdjustNorm") ; hasCustomRange = kTRUE ; } else { // Use range of last fit, if it was non-default and no other range was specified RooArgSet* plotDep = getObservables(*frame->getPlotVar()) ; RooRealVar* plotDepVar = (RooRealVar*) plotDep->find(frame->getPlotVar()->GetName()) ; if (plotDepVar->hasBinning("fit")) { rangeLo = plotDepVar->getMin("fit") ; rangeHi = plotDepVar->getMax("fit") ; adjustNorm = kTRUE ; hasCustomRange = kTRUE ; coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") INFO: pdf has been fit over restricted range, plotting only fitted " << "part of PDF normalized data in restricted range" << endl ; } delete plotDep ; } if (hasCustomRange && adjustNorm) { scaleFactor *= frame->getFitRangeNEvt(rangeLo,rangeHi)/nExpected ; } else { scaleFactor *= frame->getFitRangeNEvt()/nExpected ; } } else if (stype==RelativeExpected) { scaleFactor *= nExpected ; } else if (stype==NumEvent) { scaleFactor /= nExpected ; } scaleFactor *= frame->getFitRangeBinW() ; } frame->updateNormVars(*frame->getPlotVar()) ; // Append overriding scale factor command at end of original command list RooCmdArg tmp = RooFit::Normalization(scaleFactor,Raw) ; cmdList.Add(&tmp) ; // Was a component selected requested if (haveCompSel) { // Get complete set of tree branch nodes RooArgSet branchNodeSet ; branchNodeServerList(&branchNodeSet) ; // Discard any non-PDF nodes TIterator* iter = branchNodeSet.createIterator() ; RooAbsArg* arg ; while((arg=(RooAbsArg*)iter->Next())) { if (!dynamic_cast<RooAbsPdf*>(arg)) { branchNodeSet.remove(*arg) ; } } delete iter ; // Obtain direct selection RooArgSet* dirSelNodes ; if (compSet) { dirSelNodes = (RooArgSet*) branchNodeSet.selectCommon(*compSet) ; } else { dirSelNodes = (RooArgSet*) branchNodeSet.selectByName(compSpec) ; } coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") directly selected PDF components: " << *dirSelNodes << endl ; // Do indirect selection and activate both plotOnCompSelect(dirSelNodes) ; delete dirSelNodes ; } RooCmdArg cnsuffix("CurveNameSuffix",0,0,0,0,nameSuffix.Data(),0,0,0) ; cmdList.Add(&cnsuffix); RooPlot* ret = RooAbsReal::plotOn(frame,cmdList) ; // Restore selection status ; if (haveCompSel) plotOnCompSelect(0) ; return ret ; } //_____________________________________________________________________________ void RooAbsPdf::plotOnCompSelect(RooArgSet* selNodes) const { // Helper function for plotting of composite p.d.fs. Given // a set of selected components that should be plotted, // find all nodes that (in)directly depend on these selected // nodes. Mark all directly and indirecty selected nodes // as 'selected' using the selectComp() method // Get complete set of tree branch nodes RooArgSet branchNodeSet ; branchNodeServerList(&branchNodeSet) ; // Discard any non-PDF nodes TIterator* iter = branchNodeSet.createIterator() ; RooAbsArg* arg ; while((arg=(RooAbsArg*)iter->Next())) { if (!dynamic_cast<RooAbsPdf*>(arg)) { branchNodeSet.remove(*arg) ; } } // If no set is specified, restored all selection bits to kTRUE if (!selNodes) { // Reset PDF selection bits to kTRUE iter->Reset() ; while((arg=(RooAbsArg*)iter->Next())) { ((RooAbsPdf*)arg)->selectComp(kTRUE) ; } delete iter ; return ; } // Add all nodes below selected nodes iter->Reset() ; TIterator* sIter = selNodes->createIterator() ; RooArgSet tmp ; while((arg=(RooAbsArg*)iter->Next())) { sIter->Reset() ; RooAbsArg* selNode ; while((selNode=(RooAbsArg*)sIter->Next())) { if (selNode->dependsOn(*arg)) { tmp.add(*arg,kTRUE) ; } } } delete sIter ; // Add all nodes that depend on selected nodes iter->Reset() ; while((arg=(RooAbsArg*)iter->Next())) { if (arg->dependsOn(*selNodes)) { tmp.add(*arg,kTRUE) ; } } tmp.remove(*selNodes,kTRUE) ; tmp.remove(*this) ; selNodes->add(tmp) ; coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") indirectly selected PDF components: " << tmp << endl ; // Set PDF selection bits according to selNodes iter->Reset() ; while((arg=(RooAbsArg*)iter->Next())) { Bool_t select = selNodes->find(arg->GetName()) ? kTRUE : kFALSE ; ((RooAbsPdf*)arg)->selectComp(select) ; } delete iter ; } //_____________________________________________________________________________ RooPlot* RooAbsPdf::plotOn(RooPlot *frame, PlotOpt o) const { // Plot oneself on 'frame'. In addition to features detailed in RooAbsReal::plotOn(), // the scale factor for a PDF can be interpreted in three different ways. The interpretation // is controlled by ScaleType // // Relative - Scale factor is applied on top of PDF normalization scale factor // NumEvent - Scale factor is interpreted as a number of events. The surface area // under the PDF curve will match that of a histogram containing the specified // number of event // Raw - Scale factor is applied to the raw (projected) probability density. // Not too useful, option provided for completeness. // Sanity checks if (plotSanityChecks(frame)) return frame ; // More sanity checks Double_t nExpected(1) ; if (o.stype==RelativeExpected) { if (!canBeExtended()) { coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName() << "): ERROR the 'Expected' scale option can only be used on extendable PDFs" << endl ; return frame ; } nExpected = expectedEvents(frame->getNormVars()) ; } // Adjust normalization, if so requested if (o.stype != Raw) { if (frame->getFitRangeNEvt() && o.stype==Relative) { // If non-default plotting range is specified, adjust number of events in fit range o.scaleFactor *= frame->getFitRangeNEvt()/nExpected ; } else if (o.stype==RelativeExpected) { o.scaleFactor *= nExpected ; } else if (o.stype==NumEvent) { o.scaleFactor /= nExpected ; } o.scaleFactor *= frame->getFitRangeBinW() ; } frame->updateNormVars(*frame->getPlotVar()) ; return RooAbsReal::plotOn(frame,o) ; } //_____________________________________________________________________________ RooPlot* RooAbsPdf::paramOn(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) { // Add a box with parameter values (and errors) to the specified frame // // The following named arguments are supported // // Parameters(const RooArgSet& param) -- Only the specified subset of parameters will be shown. // By default all non-contant parameters are shown // ShowConstant(Bool_t flag) -- Also display constant parameters // Format(const char* optStr) -- Classing [arameter formatting options, provided for backward compatibility // Format(const char* what,...) -- Parameter formatting options, details given below // Label(const chat* label) -- Add header label to parameter box // Layout(Double_t xmin, -- Specify relative position of left,right side of box and top of box. Position of // Double_t xmax, Double_t ymax) bottom of box is calculated automatically from number lines in box // // // The Format(const char* what,...) has the following structure // // const char* what -- Controls what is shown. "N" adds name, "E" adds error, // "A" shows asymmetric error, "U" shows unit, "H" hides the value // FixedPrecision(int n) -- Controls precision, set fixed number of digits // AutoPrecision(int n) -- Controls precision. Number of shown digits is calculated from error // + n specified additional digits (1 is sensible default) // // Example use: pdf.paramOn(frame, Label("fit result"), Format("NEU",AutoPrecision(1)) ) ; // // Stuff all arguments in a list RooLinkedList cmdList; cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ; // Select the pdf-specific commands RooCmdConfig pc(Form("RooAbsPdf::paramOn(%s)",GetName())) ; pc.defineString("label","Label",0,"") ; pc.defineDouble("xmin","Layout",0,0.65) ; pc.defineDouble("xmax","Layout",1,0.99) ; pc.defineInt("ymaxi","Layout",0,Int_t(0.95*10000)) ; pc.defineInt("showc","ShowConstants",0,0) ; pc.defineObject("params","Parameters",0,0) ; pc.defineString("formatStr","Format",0,"NELU") ; pc.defineInt("sigDigit","Format",0,2) ; pc.defineInt("dummy","FormatArgs",0,0) ; pc.defineMutex("Format","FormatArgs") ; // Process and check varargs pc.process(cmdList) ; if (!pc.ok(kTRUE)) { return frame ; } const char* label = pc.getString("label") ; Double_t xmin = pc.getDouble("xmin") ; Double_t xmax = pc.getDouble("xmax") ; Double_t ymax = pc.getInt("ymaxi") / 10000. ; Int_t showc = pc.getInt("showc") ; const char* formatStr = pc.getString("formatStr") ; Int_t sigDigit = pc.getInt("sigDigit") ; // Decode command line arguments RooArgSet* params = static_cast<RooArgSet*>(pc.getObject("params")) ; if (!params) { params = getParameters(frame->getNormVars()) ; if (pc.hasProcessed("FormatArgs")) { const RooCmdArg* formatCmd = static_cast<RooCmdArg*>(cmdList.FindObject("FormatArgs")) ; paramOn(frame,*params,showc,label,0,0,xmin,xmax,ymax,formatCmd) ; } else { paramOn(frame,*params,showc,label,sigDigit,formatStr,xmin,xmax,ymax) ; } delete params ; } else { RooArgSet* pdfParams = getParameters(frame->getNormVars()) ; RooArgSet* selParams = static_cast<RooArgSet*>(pdfParams->selectCommon(*params)) ; if (pc.hasProcessed("FormatArgs")) { const RooCmdArg* formatCmd = static_cast<RooCmdArg*>(cmdList.FindObject("FormatArgs")) ; paramOn(frame,*selParams,showc,label,0,0,xmin,xmax,ymax,formatCmd) ; } else { paramOn(frame,*selParams,showc,label,sigDigit,formatStr,xmin,xmax,ymax) ; } delete selParams ; delete pdfParams ; } return frame ; } //_____________________________________________________________________________ RooPlot* RooAbsPdf::paramOn(RooPlot* frame, const RooAbsData* data, const char *label, Int_t sigDigits, Option_t *options, Double_t xmin, Double_t xmax ,Double_t ymax) { // OBSOLETE FUNCTION PROVIDED FOR BACKWARD COMPATIBILITY RooArgSet* params = getParameters(data) ; TString opts(options) ; paramOn(frame,*params,opts.Contains("c"),label,sigDigits,options,xmin,xmax,ymax) ; delete params ; return frame ; } //_____________________________________________________________________________ RooPlot* RooAbsPdf::paramOn(RooPlot* frame, const RooArgSet& params, Bool_t showConstants, const char *label, Int_t sigDigits, Option_t *options, Double_t xmin, Double_t xmax ,Double_t ymax, const RooCmdArg* formatCmd) { // Add a text box with the current parameter values and their errors to the frame. // Observables of this PDF appearing in the 'data' dataset will be omitted. // // Optional label will be inserted as first line of the text box. Use 'sigDigits' // to modify the default number of significant digits printed. The 'xmin,xmax,ymax' // values specify the inital relative position of the text box in the plot frame // parse the options TString opts = options; opts.ToLower(); Bool_t showLabel= (label != 0 && strlen(label) > 0); // calculate the box's size, adjusting for constant parameters TIterator* pIter = params.createIterator() ; Int_t nPar= params.getSize(); Double_t ymin(ymax), dy(0.06); Int_t index(nPar); RooRealVar *var = 0; while((var=(RooRealVar*)pIter->Next())) { if(showConstants || !var->isConstant()) ymin-= dy; } if(showLabel) ymin-= dy; // create the box and set its options TPaveText *box= new TPaveText(xmin,ymax,xmax,ymin,"BRNDC"); if(!box) return 0; box->SetFillColor(0); box->SetBorderSize(1); box->SetTextAlign(12); box->SetTextSize(0.04F); box->SetFillStyle(1001); box->SetFillColor(0); TText *text = 0; //char buffer[512]; index= nPar; pIter->Reset() ; while((var=(RooRealVar*)pIter->Next())) { if(var->isConstant() && !showConstants) continue; TString *formatted= options ? var->format(sigDigits, options) : var->format(*formatCmd) ; text= box->AddText(formatted->Data()); delete formatted; } // add the optional label if specified if(showLabel) text= box->AddText(label); // Add box to frame frame->addObject(box) ; delete pIter ; return frame ; } //_____________________________________________________________________________ Double_t RooAbsPdf::expectedEvents(const RooArgSet*) const { // Return expected number of events from this p.d.f for use in extended // likelihood calculations. This default implementation returns zero return 0 ; } //_____________________________________________________________________________ void RooAbsPdf::verboseEval(Int_t stat) { // Change global level of verbosity for p.d.f. evaluations _verboseEval = stat ; } //_____________________________________________________________________________ Int_t RooAbsPdf::verboseEval() { // Return global level of verbosity for p.d.f. evaluations return _verboseEval ; } //_____________________________________________________________________________ void RooAbsPdf::CacheElem::operModeHook(RooAbsArg::OperMode) { // Dummy implementation } //_____________________________________________________________________________ RooAbsPdf::CacheElem::~CacheElem() { // Destructor of normalization cache element. If this element // provides the 'current' normalization stored in RooAbsPdf::_norm // zero _norm pointer here before object pointed to is deleted here // Zero _norm pointer in RooAbsPdf if it is points to our cache payload if (_owner) { RooAbsPdf* pdfOwner = static_cast<RooAbsPdf*>(_owner) ; if (pdfOwner->_norm == _norm) { pdfOwner->_norm = 0 ; } } delete _norm ; } //_____________________________________________________________________________ RooAbsPdf* RooAbsPdf::createProjection(const RooArgSet& iset) { // Return a p.d.f that represent a projection of this p.d.f integrated over given observables // Construct name for new object TString name(GetName()) ; name.Append("_Proj[") ; if (iset.getSize()>0) { TIterator* iter = iset.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("]") ; // Return projected p.d.f. return new RooProjectedPdf(name.Data(),name.Data(),*this,iset) ; } //_____________________________________________________________________________ RooAbsReal* RooAbsPdf::createCdf(const RooArgSet& iset, const RooArgSet& nset) { // Create a cumulative distribution function of this p.d.f in terms // of the observables listed in iset. If no nset argument is given // the c.d.f normalization is constructed over the integrated // observables, so that its maximum value is precisely 1. It is also // possible to choose a different normalization for // multi-dimensional p.d.f.s: eg. for a pdf f(x,y,z) one can // construct a partial cdf c(x,y) that only when integrated itself // over z results in a maximum value of 1. To construct such a cdf pass // z as argument to the optional nset argument return createCdf(iset,RooFit::SupNormSet(nset)) ; } //_____________________________________________________________________________ RooAbsReal* RooAbsPdf::createCdf(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 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 // // 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::createCdf(%s)",GetName())) ; pc.defineObject("supNormSet","SupNormSet",0,0) ; pc.defineInt("numScanBins","ScanParameters",0,1000) ; pc.defineInt("intOrder","ScanParameters",1,2) ; pc.defineInt("doScanNum","ScanNumCdf",0,1) ; pc.defineInt("doScanAll","ScanAllCdf",0,0) ; pc.defineInt("doScanNon","ScanNoCdf",0,0) ; pc.defineMutex("ScanNumCdf","ScanAllCdf","ScanNoCdf") ; // 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 createScanCdf(iset,nset,numScanBins,intOrder) ; } if (doScanNum) { RooRealIntegral* tmp = (RooRealIntegral*) createIntegral(iset) ; Int_t isNum= (tmp->numIntRealVars().getSize()>0) ; delete tmp ; if (isNum) { coutI(NumIntegration) << "RooAbsPdf::createCdf(" << 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 ? createScanCdf(iset,nset,numScanBins,intOrder) : createIntRI(iset,nset) ; } return 0 ; } RooAbsReal* RooAbsPdf::createScanCdf(const RooArgSet& iset, const RooArgSet& nset, Int_t numScanBins, Int_t intOrder) { string name = string(GetName()) + "_NUMCDF_" + integralNameSuffix(iset,&nset).Data() ; RooRealVar* ivar = (RooRealVar*) iset.first() ; ivar->setBins(numScanBins,"numcdf") ; RooNumCdf* ret = new RooNumCdf(name.c_str(),name.c_str(),*this,*ivar,"numcdf") ; ret->setInterpolationOrder(intOrder) ; return ret ; } //_____________________________________________________________________________ RooArgSet* RooAbsPdf::getAllConstraints(const RooArgSet& observables, const RooArgSet& constrainedParams) const { // This helper function finds and collects all constraints terms of all coponent p.d.f.s // and returns a RooArgSet with all those terms RooArgSet* ret = new RooArgSet("AllConstraints") ; RooArgSet* comps = getComponents() ; TIterator* iter = comps->createIterator() ; RooAbsArg *arg ; while((arg=(RooAbsArg*)iter->Next())) { RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ; if (pdf) { RooArgSet* compRet = pdf->getConstraints(observables,constrainedParams) ; if (compRet) { ret->add(*compRet,kFALSE) ; delete compRet ; } } } delete iter ; delete comps ; return ret ; } //_____________________________________________________________________________ void RooAbsPdf::clearEvalError() { // Clear the evaluation error flag _evalError = kFALSE ; } //_____________________________________________________________________________ Bool_t RooAbsPdf::evalError() { // Return the evaluation error flag return _evalError ; } //_____________________________________________________________________________ Bool_t RooAbsPdf::isSelectedComp() const { // If true, the current pdf is a selected component (for use in plotting) return _selectComp || _globalSelectComp ; } //_____________________________________________________________________________ void RooAbsPdf::globalSelectComp(Bool_t flag) { // Global switch controlling the activation of the selectComp() functionality _globalSelectComp = flag ; } //_____________________________________________________________________________ void RooAbsPdf::raiseEvalError() { // Raise the evaluation error flag _evalError = kTRUE ; }