Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooAbsReal.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17//////////////////////////////////////////////////////////////////////////////
18
19/** \class RooAbsReal
20
21 Abstract base class for objects that represent a
22 real value and implements functionality common to all real-valued objects
23 such as the ability to plot them, to construct integrals of them, the
24 ability to advertise (partial) analytical integrals etc.
25
26 Implementation of RooAbsReal may be derived, thus no interface
27 is provided to modify the contents.
28
29 \ingroup Roofitcore
30*/
31
32#include "RooAbsReal.h"
33
34#include "FitHelpers.h"
36#include "RooAbsData.h"
37#include "RooAddPdf.h"
38#include "RooAddition.h"
39#include "RooArgList.h"
40#include "RooArgProxy.h"
41#include "RooArgSet.h"
42#include "RooBinning.h"
43#include "RooBrentRootFinder.h"
44#include "RooCachedReal.h"
45#include "RooCategory.h"
46#include "RooCmdConfig.h"
47#include "RooConstVar.h"
48#include "RooCurve.h"
49#include "RooCustomizer.h"
50#include "RooDataHist.h"
51#include "RooDataSet.h"
52#include "RooDerivative.h"
53#include "RooFirstMoment.h"
55#include "RooFit/Evaluator.h"
56#include "RooFitResult.h"
57#include "RooFormulaVar.h"
58#include "RooFunctor.h"
59#include "RooGlobalFunc.h"
60#include "RooFitImplHelpers.h"
61#include "RooHist.h"
62#include "RooMoment.h"
63#include "RooMsgService.h"
64#include "RooNumIntConfig.h"
65#include "RooNumRunningInt.h"
66#include "RooParamBinning.h"
67#include "RooPlot.h"
68#include "RooProduct.h"
69#include "RooProfileLL.h"
70#include "RooRealBinding.h"
71#include "RooRealIntegral.h"
72#include "RooRealVar.h"
73#include "RooSecondMoment.h"
74#include "RooVectorDataStore.h"
75#include "TreeReadBuffer.h"
76#include "ValueChecking.h"
77
78#include "ROOT/StringUtils.hxx"
79#include "Compression.h"
80#include "Math/IFunction.h"
81#include "TMath.h"
82#include "TObjString.h"
83#include "TTree.h"
84#include "TH1.h"
85#include "TH2.h"
86#include "TH3.h"
87#include "TBranch.h"
88#include "TLeaf.h"
89#include "TAttLine.h"
90#include "TF1.h"
91#include "TF2.h"
92#include "TF3.h"
93#include "TMatrixD.h"
94#include "TVector.h"
95#include "strlcpy.h"
96#ifndef NDEBUG
97#include <TSystem.h> // To print stack traces when caching errors are detected
98#endif
99
100#include <iomanip>
101#include <iostream>
102#include <limits>
103#include <sstream>
104#include <sys/types.h>
105
106namespace {
107
108// Internal helper RooAbsFunc that evaluates the scaled data-weighted average of
109// given RooAbsReal as a function of a single variable using the RooFit::Evaluator.
110class ScaledDataWeightedAverage : public RooAbsFunc {
111public:
112 ScaledDataWeightedAverage(RooAbsReal const &arg, RooAbsData const &data, double scaleFactor, RooAbsRealLValue &var)
113 : RooAbsFunc{1}, _var{var}, _dataWeights{data.getWeightBatch(0, data.numEntries())}, _scaleFactor{scaleFactor}
114 {
115 _arg = RooFit::Detail::compileForNormSet(arg, *data.get());
116 _arg->recursiveRedirectServers(RooArgList{var});
117 _evaluator = std::make_unique<RooFit::Evaluator>(*_arg);
118 std::stack<std::vector<double>>{}.swap(_vectorBuffers);
119 auto dataSpans =
120 RooFit::Detail::BatchModeDataHelpers::getDataSpans(data, "", nullptr, /*skipZeroWeights=*/false,
121 /*takeGlobalObservablesFromData=*/true, _vectorBuffers);
122 for (auto const& item : dataSpans) {
123 _evaluator->setInput(item.first->GetName(), item.second, false);
124 }
125 }
126
127 double operator()(const double xvector[]) const override
128 {
129 double oldVal = _var.getVal();
130 _var.setVal(xvector[0]);
131
132 double out = 0.0;
133 std::span<const double> pdfValues = _evaluator->run();
134 if (_dataWeights.empty()) {
135 out = std::accumulate(pdfValues.begin(), pdfValues.end(), 0.0) / pdfValues.size();
136 } else {
137 double weightsSum = 0.0;
138 for (std::size_t i = 0; i < pdfValues.size(); ++i) {
139 out += pdfValues[i] * _dataWeights[i];
140 weightsSum += _dataWeights[i];
141 }
142 out /= weightsSum;
143 }
144 out *= _scaleFactor;
145
146 _var.setVal(oldVal);
147 return out;
148 }
149 double getMinLimit(UInt_t /*dimension*/) const override { return _var.getMin(); }
150 double getMaxLimit(UInt_t /*dimension*/) const override { return _var.getMax(); }
151
152private:
153 RooAbsRealLValue &_var;
154 std::unique_ptr<RooAbsReal> _arg;
155 std::span<const double> _dataWeights;
156 double _scaleFactor;
157 std::unique_ptr<RooFit::Evaluator> _evaluator;
158 std::stack<std::vector<double>> _vectorBuffers;
159};
160
161} // namespace
162
163using namespace std ;
164
166
168bool RooAbsReal::_hideOffset = true ;
169
170void RooAbsReal::setHideOffset(bool flag) { _hideOffset = flag ; }
172
175std::map<const RooAbsArg*,std::pair<std::string,std::list<RooAbsReal::EvalError> > > RooAbsReal::_evalErrorList ;
176
177
178////////////////////////////////////////////////////////////////////////////////
179/// coverity[UNINIT_CTOR]
180/// Default constructor
181
183
184
185////////////////////////////////////////////////////////////////////////////////
186/// Constructor with unit label
187
188RooAbsReal::RooAbsReal(const char *name, const char *title, const char *unit) :
189 RooAbsArg(name,title), _unit(unit)
190{
191 setValueDirty() ;
192 setShapeDirty() ;
193}
194
195
196////////////////////////////////////////////////////////////////////////////////
197/// Constructor with plot range and unit label
198
199RooAbsReal::RooAbsReal(const char *name, const char *title, double inMinVal,
200 double inMaxVal, const char *unit) :
201 RooAbsArg(name,title), _plotMin(inMinVal), _plotMax(inMaxVal), _unit(unit)
202{
203 setValueDirty() ;
204 setShapeDirty() ;
205}
206
207
208////////////////////////////////////////////////////////////////////////////////
209/// Copy constructor
210RooAbsReal::RooAbsReal(const RooAbsReal& other, const char* name) :
211 RooAbsArg(other,name), _plotMin(other._plotMin), _plotMax(other._plotMax),
212 _plotBins(other._plotBins), _value(other._value), _unit(other._unit), _label(other._label),
213 _forceNumInt(other._forceNumInt), _selectComp(other._selectComp)
214{
215 if (other._specIntegratorConfig) {
216 _specIntegratorConfig = std::make_unique<RooNumIntConfig>(*other._specIntegratorConfig) ;
217 }
218}
219
220
221////////////////////////////////////////////////////////////////////////////////
222/// Destructor
223
225
226
227////////////////////////////////////////////////////////////////////////////////
228/// Equality operator comparing to a double
229
231{
232 return (getVal()==value) ;
233}
234
235
236
237////////////////////////////////////////////////////////////////////////////////
238/// Equality operator when comparing to another RooAbsArg.
239/// Only functional when the other arg is a RooAbsReal
240
241bool RooAbsReal::operator==(const RooAbsArg& other) const
242{
243 const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
244 return otherReal ? operator==(otherReal->getVal()) : false ;
245}
246
247
248////////////////////////////////////////////////////////////////////////////////
249
250bool RooAbsReal::isIdentical(const RooAbsArg& other, bool assumeSameType) const
251{
252 if (!assumeSameType) {
253 const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
254 return otherReal ? operator==(otherReal->getVal()) : false ;
255 } else {
256 return getVal() == static_cast<const RooAbsReal&>(other).getVal();
257 }
258}
259
260
261////////////////////////////////////////////////////////////////////////////////
262/// Return this variable's title string. If appendUnit is true and
263/// this variable has units, also append a string " (<unit>)".
264
265TString RooAbsReal::getTitle(bool appendUnit) const
266{
267 TString title(GetTitle());
268 if(appendUnit && 0 != strlen(getUnit())) {
269 title.Append(" (");
270 title.Append(getUnit());
271 title.Append(")");
272 }
273 return title;
274}
275
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// Return value of object. If the cache is clean, return the
280/// cached value, otherwise recalculate on the fly and refill
281/// the cache
282
283double RooAbsReal::getValV(const RooArgSet* nset) const
284{
285 if (nset && nset->uniqueId().value() != _lastNormSetId) {
286 const_cast<RooAbsReal*>(this)->setProxyNormSet(nset);
287 _lastNormSetId = nset->uniqueId().value();
288 }
289
290 if (isValueDirtyAndClear()) {
291 _value = traceEval(nullptr) ;
292 // clearValueDirty() ;
293 }
294 // cout << "RooAbsReal::getValV(" << GetName() << ") writing _value = " << _value << std::endl ;
295
296 return hideOffset() ? _value + offset() : _value;
297}
298
299
300////////////////////////////////////////////////////////////////////////////////
301
303{
304 return _evalErrorList.size() ;
305}
306
307
308////////////////////////////////////////////////////////////////////////////////
309
311{
312 return _evalErrorList.begin() ;
313}
314
315
316////////////////////////////////////////////////////////////////////////////////
317/// Calculate current value of object, with error tracing wrapper
318
319double RooAbsReal::traceEval(const RooArgSet* /*nset*/) const
320{
321 double value = evaluate() ;
322
323 if (TMath::IsNaN(value)) {
324 logEvalError("function value is NAN") ;
325 }
326
327 //cxcoutD(Tracing) << "RooAbsReal::getValF(" << GetName() << ") operMode = " << _operMode << " recalculated, new value = " << value << std::endl ;
328
329 //Standard tracing code goes here
330 if (!isValidReal(value)) {
331 coutW(Tracing) << "RooAbsReal::traceEval(" << GetName()
332 << "): validation failed: " << value << std::endl ;
333 }
334
335 //Call optional subclass tracing code
336 // traceEvalHook(value) ;
337
338 return value ;
339}
340
341
342
343////////////////////////////////////////////////////////////////////////////////
344/// Variant of getAnalyticalIntegral that is also passed the normalization set
345/// that should be applied to the integrand of which the integral is requested.
346/// For certain operator p.d.f it is useful to overload this function rather
347/// than analyticalIntegralWN() as the additional normalization information
348/// may be useful in determining a more efficient decomposition of the
349/// requested integral.
350
352 const RooArgSet* /*normSet*/, const char* rangeName) const
353{
354 return _forceNumInt ? 0 : getAnalyticalIntegral(allDeps,analDeps,rangeName) ;
355}
356
357
358
359////////////////////////////////////////////////////////////////////////////////
360/// Interface function getAnalyticalIntergral advertises the
361/// analytical integrals that are supported. 'integSet'
362/// is the set of dependents for which integration is requested. The
363/// function should copy the subset of dependents it can analytically
364/// integrate to anaIntSet and return a unique identification code for
365/// this integration configuration. If no integration can be
366/// performed, zero should be returned.
367
368Int_t RooAbsReal::getAnalyticalIntegral(RooArgSet& /*integSet*/, RooArgSet& /*anaIntSet*/, const char* /*rangeName*/) const
369{
370 return 0 ;
371}
372
373
374
375////////////////////////////////////////////////////////////////////////////////
376/// Implements the actual analytical integral(s) advertised by
377/// getAnalyticalIntegral. This functions will only be called with
378/// codes returned by getAnalyticalIntegral, except code zero.
379
380double RooAbsReal::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
381{
382// cout << "RooAbsReal::analyticalIntegralWN(" << GetName() << ") code = " << code << " normSet = " << (normSet?*normSet:RooArgSet()) << std::endl ;
383 if (code==0) return getVal(normSet) ;
384 return analyticalIntegral(code,rangeName) ;
385}
386
387
388
389////////////////////////////////////////////////////////////////////////////////
390/// Implements the actual analytical integral(s) advertised by
391/// getAnalyticalIntegral. This functions will only be called with
392/// codes returned by getAnalyticalIntegral, except code zero.
393
394double RooAbsReal::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
395{
396 // By default no analytical integrals are implemented
397 coutF(Eval) << "RooAbsReal::analyticalIntegral(" << GetName() << ") code " << code << " not implemented" << std::endl ;
398 return 0 ;
399}
400
401
402
403////////////////////////////////////////////////////////////////////////////////
404/// Get the label associated with the variable
405
406const char *RooAbsReal::getPlotLabel() const
407{
408 return _label.IsNull() ? fName.Data() : _label.Data();
409}
410
411
412
413////////////////////////////////////////////////////////////////////////////////
414/// Set the label associated with this variable
415
416void RooAbsReal::setPlotLabel(const char *label)
417{
418 _label= label;
419}
420
421
422
423////////////////////////////////////////////////////////////////////////////////
424///Read object contents from stream (dummy for now)
425
426bool RooAbsReal::readFromStream(std::istream& /*is*/, bool /*compact*/, bool /*verbose*/)
427{
428 return false ;
429}
430
431
432
433////////////////////////////////////////////////////////////////////////////////
434///Write object contents to stream (dummy for now)
435
436void RooAbsReal::writeToStream(std::ostream& /*os*/, bool /*compact*/) const
437{
438}
439
440
441
442////////////////////////////////////////////////////////////////////////////////
443/// Print object value
444
445void RooAbsReal::printValue(std::ostream& os) const
446{
447 os << getVal() ;
448}
449
450
451
452////////////////////////////////////////////////////////////////////////////////
453/// Structure printing
454
455void RooAbsReal::printMultiline(std::ostream& os, Int_t contents, bool verbose, TString indent) const
456{
457 RooAbsArg::printMultiline(os,contents,verbose,indent) ;
458 os << indent << "--- RooAbsReal ---" << std::endl;
459 TString unit(_unit);
460 if(!unit.IsNull()) unit.Prepend(' ');
461 //os << indent << " Value = " << getVal() << unit << std::endl;
462 os << std::endl << indent << " Plot label is \"" << getPlotLabel() << "\"" << "\n";
463}
464
465
466////////////////////////////////////////////////////////////////////////////////
467/// Create a RooProfileLL object that eliminates all nuisance parameters in the
468/// present function. The nuisance parameters are defined as all parameters
469/// of the function except the stated paramsOfInterest
470
472{
473 // Construct name of profile object
474 auto name = std::string(GetName()) + "_Profile[";
475 bool first = true;
476 for (auto const& arg : paramsOfInterest) {
477 if (first) {
478 first = false ;
479 } else {
480 name.append(",") ;
481 }
482 name.append(arg->GetName()) ;
483 }
484 name.append("]") ;
485
486 // Create and return profile object
487 auto out = std::make_unique<RooProfileLL>(name.c_str(),(std::string("Profile of ") + GetTitle()).c_str(),*this,paramsOfInterest);
488 return RooFit::makeOwningPtr<RooAbsReal>(std::move(out));
489}
490
491
492
493
494
495
496////////////////////////////////////////////////////////////////////////////////
497/// Create an object that represents the integral of the function over one or more observables listed in `iset`.
498/// The actual integration calculation is only performed when the returned object is evaluated. The name
499/// of the integral object is automatically constructed from the name of the input function, the variables
500/// it integrates and the range integrates over.
501///
502/// \note The integral over a PDF is usually not normalised (*i.e.*, it is usually not
503/// 1 when integrating the PDF over the full range). In fact, this integral is used *to compute*
504/// the normalisation of each PDF. See the [rf110 tutorial](group__tutorial__roofit.html)
505/// for details on PDF normalisation.
506///
507/// The following named arguments are accepted
508/// | | Effect on integral creation
509/// |--|-------------------------------
510/// | `NormSet(const RooArgSet&)` | Specify normalization set, mostly useful when working with PDFs
511/// | `NumIntConfig(const RooNumIntConfig&)` | Use given configuration for any numeric integration, if necessary
512/// | `Range(const char* name)` | Integrate only over given range. Multiple ranges may be specified by passing multiple Range() arguments
513
515 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
516 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
517{
518 // Define configuration for this method
519 RooCmdConfig pc("RooAbsReal::createIntegral(" + std::string(GetName()) + ")");
520 pc.defineString("rangeName","RangeWithName",0,"",true) ;
521 pc.defineSet("normSet","NormSet",0,nullptr) ;
522 pc.defineObject("numIntConfig","NumIntConfig",0,nullptr) ;
523
524 // Process & check varargs
525 pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
526 if (!pc.ok(true)) {
527 return nullptr;
528 }
529
530 // Extract values from named arguments
531 const char* rangeName = pc.getString("rangeName",nullptr,true) ;
532 const RooArgSet* nset = pc.getSet("normSet",nullptr);
533 const RooNumIntConfig* cfg = static_cast<const RooNumIntConfig*>(pc.getObject("numIntConfig",nullptr)) ;
534
535 return createIntegral(iset,nset,cfg,rangeName) ;
536}
537
538
539
540
541
542////////////////////////////////////////////////////////////////////////////////
543/// Create an object that represents the integral of the function over one or more observables listed in iset.
544/// The actual integration calculation is only performed when the return object is evaluated. The name
545/// of the integral object is automatically constructed from the name of the input function, the variables
546/// it integrates and the range integrates over. If nset is specified the integrand is request
547/// to be normalized over nset (only meaningful when the integrand is a pdf). If rangename is specified
548/// the integral is performed over the named range, otherwise it is performed over the domain of each
549/// integrated observable. If cfg is specified it will be used to configure any numeric integration
550/// aspect of the integral. It will not force the integral to be performed numerically, which is
551/// decided automatically by RooRealIntegral.
552
554 const RooNumIntConfig* cfg, const char* rangeName) const
555{
556 if (!rangeName || strchr(rangeName,',')==nullptr) {
557 // Simple case: integral over full range or single limited range
558 return createIntObj(iset,nset,cfg,rangeName);
559 }
560
561 // Integral over multiple ranges
562 std::vector<std::string> tokens = ROOT::Split(rangeName, ",");
563
564 if(RooHelpers::checkIfRangesOverlap(iset, tokens)) {
565 std::stringstream errMsg;
566 errMsg << GetName() << " : integrating with respect to the variables " << iset << " on the ranges \"" << rangeName
567 << "\" is not possible because the ranges are overlapping";
568 const std::string errMsgString = errMsg.str();
569 coutE(Integration) << errMsgString << std::endl;
570 throw std::invalid_argument(errMsgString);
571 }
572
573 RooArgSet components ;
574 for (const std::string& token : tokens) {
575 components.addOwned(std::unique_ptr<RooAbsReal>{createIntObj(iset,nset,cfg, token.c_str())});
576 }
577
578 const std::string title = std::string("Integral of ") + GetTitle();
579 const std::string fullName = std::string(GetName()) + integralNameSuffix(iset,nset,rangeName).Data();
580
581 auto out = std::make_unique<RooAddition>(fullName.c_str(), title.c_str(), components);
582 out->addOwnedComponents(std::move(components));
583 return RooFit::makeOwningPtr<RooAbsReal>(std::move(out));
584}
585
586
587
588////////////////////////////////////////////////////////////////////////////////
589/// Internal utility function for createIntegral() that creates the actual integral object.
591 const RooNumIntConfig* cfg, const char* rangeName) const
592{
593 // Make internal use copies of iset and nset
594 RooArgSet iset(iset2) ;
595 const RooArgSet* nset = nset2 ;
596
597
598 // Initialize local variables perparing for recursive loop
599 bool error = false ;
600 const RooAbsReal* integrand = this ;
601 std::unique_ptr<RooAbsReal> integral;
602
603 // Handle trivial case of no integration here explicitly
604 if (iset.empty()) {
605
606 const std::string title = std::string("Integral of ") + GetTitle();
607 const std::string name = std::string(GetName()) + integralNameSuffix(iset,nset,rangeName).Data();
608
609 auto out = std::make_unique<RooRealIntegral>(name.c_str(), title.c_str(), *this, iset, nset, cfg, rangeName);
610 return RooFit::makeOwningPtr<RooAbsReal>(std::move(out));
611 }
612
613 // Process integration over remaining integration variables
614 while(!iset.empty()) {
615
616
617 // Find largest set of observables that can be integrated in one go
618 RooArgSet innerSet ;
619 findInnerMostIntegration(iset,innerSet,rangeName) ;
620
621 // If largest set of observables that can be integrated is empty set, problem was ill defined
622 // Postpone error messaging and handling to end of function, exit loop here
623 if (innerSet.empty()) {
624 error = true ;
625 break ;
626 }
627
628 // Prepare name and title of integral to be created
629 const std::string title = std::string("Integral of ") + integrand->GetTitle();
630 const std::string name = std::string(integrand->GetName()) + integrand->integralNameSuffix(innerSet,nset,rangeName).Data();
631
632 std::unique_ptr<RooAbsReal> innerIntegral = std::move(integral);
633
634 // Construct innermost integral
635 integral = std::make_unique<RooRealIntegral>(name.c_str(),title.c_str(),*integrand,innerSet,nset,cfg,rangeName);
636
637 // Integral of integral takes ownership of innermost integral
638 if (innerIntegral) {
639 integral->addOwnedComponents(std::move(innerIntegral));
640 }
641
642 // Remove already integrated observables from to-do list
643 iset.remove(innerSet) ;
644
645 // Send info message on recursion if needed
646 if (integrand == this && !iset.empty()) {
647 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" << std::endl ;
648 }
649
650 // Prepare for recursion, next integral should integrate last integrand
651 integrand = integral.get();
652
653
654 // Only need normalization set in innermost integration
655 nset = nullptr;
656 }
657
658 if (error) {
659 coutE(Integration) << GetName() << " : ERROR while defining recursive integral over observables with parameterized integration ranges, please check that integration rangs specify uniquely defined integral " << std::endl;
660 return nullptr;
661 }
662
663
664 // After-burner: apply interpolating cache on (numeric) integral if requested by user
665 const char* cacheParamsStr = getStringAttribute("CACHEPARAMINT") ;
666 if (cacheParamsStr && strlen(cacheParamsStr)) {
667
668 std::unique_ptr<RooArgSet> intParams{integral->getVariables()};
669
670 RooArgSet cacheParams = RooHelpers::selectFromArgSet(*intParams, cacheParamsStr);
671
672 if (!cacheParams.empty()) {
673 cxcoutD(Caching) << "RooAbsReal::createIntObj(" << GetName() << ") INFO: constructing " << cacheParams.size()
674 << "-dim value cache for integral over " << iset2 << " as a function of " << cacheParams << " in range " << (rangeName?rangeName:"<none>") << std::endl ;
675 std::string name = Form("%s_CACHE_[%s]",integral->GetName(),cacheParams.contentsString().c_str()) ;
676 auto cachedIntegral = std::make_unique<RooCachedReal>(name.c_str(),name.c_str(),*integral,cacheParams);
677 cachedIntegral->setInterpolationOrder(2) ;
678 cachedIntegral->addOwnedComponents(std::move(integral));
679 cachedIntegral->setCacheSource(true) ;
680 if (integral->operMode()==ADirty) {
681 cachedIntegral->setOperMode(ADirty) ;
682 }
683 //cachedIntegral->disableCache(true) ;
684 return RooFit::makeOwningPtr<RooAbsReal>(std::move(cachedIntegral));
685 }
686 }
687
688 return RooFit::makeOwningPtr(std::move(integral));
689}
690
691
692
693////////////////////////////////////////////////////////////////////////////////
694/// Utility function for createIntObj() that aids in the construct of recursive integrals
695/// over functions with multiple observables with parameterized ranges. This function
696/// finds in a given set allObs over which integration is requested the largeset subset
697/// of observables that can be integrated simultaneously. This subset consists of
698/// observables with fixed ranges and observables with parameterized ranges whose
699/// parameterization does not depend on any observable that is also integrated.
700
701void RooAbsReal::findInnerMostIntegration(const RooArgSet& allObs, RooArgSet& innerObs, const char* rangeName) const
702{
703 // Make lists of
704 // a) integrated observables with fixed ranges,
705 // b) integrated observables with parameterized ranges depending on other integrated observables
706 // c) integrated observables used in definition of any parameterized ranges of integrated observables
707 RooArgSet obsWithFixedRange(allObs) ;
708 RooArgSet obsWithParamRange ;
709 RooArgSet obsServingAsRangeParams ;
710
711 // Loop over all integrated observables
712 for (const auto aarg : allObs) {
713 // Check if observable is real-valued lvalue
714 if (auto arglv = dynamic_cast<RooAbsRealLValue*>(aarg)) {
715
716 // Check if range is parameterized
717 RooAbsBinning& binning = arglv->getBinning(rangeName,false,true) ;
718 if (binning.isParameterized()) {
719 RooArgSet loBoundObs;
720 RooArgSet hiBoundObs;
721 binning.lowBoundFunc()->getObservables(&allObs, loBoundObs) ;
722 binning.highBoundFunc()->getObservables(&allObs, hiBoundObs) ;
723
724 // Check if range parameterization depends on other integrated observables
725 if (loBoundObs.overlaps(allObs) || hiBoundObs.overlaps(allObs)) {
726 obsWithParamRange.add(*aarg) ;
727 obsWithFixedRange.remove(*aarg) ;
728 obsServingAsRangeParams.add(loBoundObs,false) ;
729 obsServingAsRangeParams.add(hiBoundObs,false) ;
730 }
731 }
732 }
733 }
734
735 // Make list of fixed-range observables that are _not_ involved in the parameterization of ranges of other observables
736 RooArgSet obsWithFixedRangeNP(obsWithFixedRange) ;
737 obsWithFixedRangeNP.remove(obsServingAsRangeParams) ;
738
739 // Make list of param-range observables that are _not_ involved in the parameterization of ranges of other observables
740 RooArgSet obsWithParamRangeNP(obsWithParamRange) ;
741 obsWithParamRangeNP.remove(obsServingAsRangeParams) ;
742
743 // Construct inner-most integration: over observables (with fixed or param range) not used in any other param range definitions
744 innerObs.removeAll() ;
745 innerObs.add(obsWithFixedRangeNP) ;
746 innerObs.add(obsWithParamRangeNP) ;
747
748}
749
750
751////////////////////////////////////////////////////////////////////////////////
752/// Construct string with unique suffix name to give to integral object that encodes
753/// integrated observables, normalization observables and the integration range name
754
755TString RooAbsReal::integralNameSuffix(const RooArgSet& iset, const RooArgSet* nset, const char* rangeName, bool omitEmpty) const
756{
757 TString name ;
758 if (!iset.empty()) {
759
760 RooArgSet isetTmp(iset) ;
761 isetTmp.sort() ;
762
763 name.Append("_Int[") ;
764 bool first(true) ;
765 for(RooAbsArg * arg : isetTmp) {
766 if (first) {
767 first=false ;
768 } else {
769 name.Append(",") ;
770 }
771 name.Append(arg->GetName()) ;
772 }
773 if (rangeName) {
774 name.Append("|") ;
775 name.Append(rangeName) ;
776 }
777 name.Append("]");
778 } else if (!omitEmpty) {
779 name.Append("_Int[]") ;
780 }
781
782 if (nset && !nset->empty()) {
783
784 RooArgSet nsetTmp(*nset) ;
785 nsetTmp.sort() ;
786
787 name.Append("_Norm[") ;
788 bool first(true);
789 for(RooAbsArg * arg : nsetTmp) {
790 if (first) {
791 first=false ;
792 } else {
793 name.Append(",") ;
794 }
795 name.Append(arg->GetName()) ;
796 }
797 const RooAbsPdf* thisPdf = dynamic_cast<const RooAbsPdf*>(this) ;
798 if (thisPdf && thisPdf->normRange()) {
799 name.Append("|") ;
800 name.Append(thisPdf->normRange()) ;
801 }
802 name.Append("]") ;
803 }
804
805 return name ;
806}
807
808
809
810////////////////////////////////////////////////////////////////////////////////
811/// Utility function for plotOn() that creates a projection of a function or p.d.f
812/// to be plotted on a RooPlot.
813/// \ref createPlotProjAnchor "createPlotProjection()"
814
815const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars,
816 RooArgSet*& cloneSet) const
817{
818 return createPlotProjection(depVars,&projVars,cloneSet) ;
819}
820
821
822////////////////////////////////////////////////////////////////////////////////
823/// Utility function for plotOn() that creates a projection of a function or p.d.f
824/// to be plotted on a RooPlot.
825/// \anchor createPlotProjAnchor
826///
827/// Create a new object \f$ G \f$ that represents the normalized projection:
828/// \f[
829/// G[x,p] = \frac{\int F[x,y,p] \; \mathrm{d}\{y\}}
830/// {\int F[x,y,p] \; \mathrm{d}\{x\} \, \mathrm{d}\{y\}}
831/// \f]
832/// where \f$ F[x,y,p] \f$ is the function we represent, and
833/// \f$ \{ p \} \f$ are the remaining variables ("parameters").
834///
835/// \param[in] dependentVars Dependent variables over which to normalise, \f$ \{x\} \f$.
836/// \param[in] projectedVars Variables to project out, \f$ \{ y \} \f$.
837/// \param[out] cloneSet Will be set to a RooArgSet*, which will contain a clone of *this plus its projection integral object.
838/// The latter will also be returned. The caller takes ownership of this set.
839/// \param[in] rangeName Optional range for projection integrals
840/// \param[in] condObs Conditional observables, which are not integrated for normalisation, even if they
841/// are in `dependentVars` or `projectedVars`.
842/// \return A pointer to the newly created object, or zero in case of an
843/// error. The caller is responsible for deleting the `cloneSet` (which includes the returned projection object).
844const RooAbsReal *RooAbsReal::createPlotProjection(const RooArgSet &dependentVars, const RooArgSet *projectedVars,
845 RooArgSet *&cloneSet, const char* rangeName, const RooArgSet* condObs) const
846{
847 // Get the set of our leaf nodes
848 RooArgSet leafNodes;
849 RooArgSet treeNodes;
850 leafNodeServerList(&leafNodes,this);
851 treeNodeServerList(&treeNodes,this) ;
852
853
854 // Check that the dependents are all fundamental. Filter out any that we
855 // do not depend on, and make substitutions by name in our leaf list.
856 // Check for overlaps with the projection variables.
857 for (const auto arg : dependentVars) {
858 if(!arg->isFundamental() && !dynamic_cast<const RooAbsLValue*>(arg)) {
859 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: variable \"" << arg->GetName()
860 << "\" of wrong type: " << arg->ClassName() << std::endl;
861 return nullptr;
862 }
863
864 RooAbsArg *found= treeNodes.find(arg->GetName());
865 if(!found) {
866 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
867 << "\" is not a dependent and will be ignored." << std::endl;
868 continue;
869 }
870 if(found != arg) {
871 if (leafNodes.find(found->GetName())) {
872 leafNodes.replace(*found,*arg);
873 } else {
874 leafNodes.add(*arg) ;
875
876 // Remove any dependents of found, replace by dependents of LV node
877 RooArgSet lvDep;
878 arg->getObservables(&leafNodes, lvDep);
879 for (const auto lvs : lvDep) {
880 RooAbsArg* tmp = leafNodes.find(lvs->GetName()) ;
881 if (tmp) {
882 leafNodes.remove(*tmp) ;
883 leafNodes.add(*lvs) ;
884 }
885 }
886 }
887 }
888
889 // check if this arg is also in the projection set
890 if(nullptr != projectedVars && projectedVars->find(arg->GetName())) {
891 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
892 << "\" cannot be both a dependent and a projected variable." << std::endl;
893 return nullptr;
894 }
895 }
896
897 // Remove the projected variables from the list of leaf nodes, if necessary.
898 if(nullptr != projectedVars) leafNodes.remove(*projectedVars,true);
899
900 // Make a deep-clone of ourself so later operations do not disturb our original state
901 cloneSet = new RooArgSet;
902 if (RooArgSet(*this).snapshot(*cloneSet, true)) {
903 coutE(Plotting) << "RooAbsPdf::createPlotProjection(" << GetName() << ") Couldn't deep-clone PDF, abort," << std::endl ;
904 return nullptr ;
905 }
906 RooAbsReal *theClone= static_cast<RooAbsReal*>(cloneSet->find(GetName()));
907
908 // The remaining entries in our list of leaf nodes are the external
909 // dependents (x) and parameters (p) of the projection. Patch them back
910 // into the theClone. This orphans the nodes they replace, but the orphans
911 // are still in the cloneList and so will be cleaned up eventually.
912 //cout << "redirection leafNodes : " ; leafNodes.Print("1") ;
913
914 std::unique_ptr<RooArgSet> plotLeafNodes{static_cast<RooArgSet*>(leafNodes.selectCommon(dependentVars))};
915 theClone->recursiveRedirectServers(*plotLeafNodes,false,false,false);
916
917 // Create the set of normalization variables to use in the projection integrand
918 RooArgSet normSet(dependentVars);
919 if(nullptr != projectedVars) normSet.add(*projectedVars);
920 if(nullptr != condObs) {
921 normSet.remove(*condObs,true,true) ;
922 }
923
924 // Try to create a valid projection integral. If no variables are to be projected,
925 // create a null projection anyway to bind our normalization over the dependents
926 // consistently with the way they would be bound with a non-trivial projection.
927 RooArgSet empty;
928 if(nullptr == projectedVars) projectedVars= &empty;
929
930 TString name = GetName() ;
931 name += integralNameSuffix(*projectedVars,&normSet,rangeName,true) ;
932
933 TString title(GetTitle());
934 title.Prepend("Projection of ");
935
936
937 std::unique_ptr<RooAbsReal> projected{theClone->createIntegral(*projectedVars,normSet,rangeName)};
938
939 if(nullptr == projected || !projected->isValid()) {
940 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: cannot integrate out ";
941 projectedVars->printStream(std::cout,kName|kArgs,kSingleLine);
942 return nullptr;
943 }
944
945 if(projected->InheritsFrom(RooRealIntegral::Class())){
946 static_cast<RooRealIntegral&>(*projected).setAllowComponentSelection(true);
947 }
948
949 projected->SetName(name.Data()) ;
950 projected->SetTitle(title.Data()) ;
951
952 // Add the projection integral to the cloneSet so that it eventually gets cleaned up by the caller.
953 RooAbsReal *projectedPtr = projected.get();
954 cloneSet->addOwned(std::move(projected));
955
956 // return a const pointer to remind the caller that they do not delete the returned object
957 // directly (it is contained in the cloneSet instead).
958 return projectedPtr;
959}
960
961
962
963////////////////////////////////////////////////////////////////////////////////
964/// Fill the ROOT histogram 'hist' with values sampled from this
965/// function at the bin centers. Our value is calculated by first
966/// integrating out any variables in projectedVars and then scaling
967/// the result by scaleFactor. Returns a pointer to the input
968/// histogram, or zero in case of an error. The input histogram can
969/// be any TH1 subclass, and therefore of arbitrary
970/// dimension. Variables are matched with the (x,y,...) dimensions of
971/// the input histogram according to the order in which they appear
972/// in the input plotVars list. If scaleForDensity is true the
973/// histogram is filled with a the functions density rather than
974/// the functions value (i.e. the value at the bin center is multiplied
975/// with bin volume)
976
978 double scaleFactor, const RooArgSet *projectedVars, bool scaleForDensity,
979 const RooArgSet* condObs, bool setError) const
980{
981 // Do we have a valid histogram to use?
982 if(nullptr == hist) {
983 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: no valid histogram to fill" << std::endl;
984 return nullptr;
985 }
986
987 // Check that the number of plotVars matches the input histogram's dimension
988 Int_t hdim= hist->GetDimension();
989 if(hdim != int(plotVars.size())) {
990 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: plotVars has the wrong dimension" << std::endl;
991 return nullptr;
992 }
993
994
995 // Check that the plot variables are all actually RooRealVars and print a warning if we do not
996 // explicitly depend on one of them. Fill a set (not list!) of cloned plot variables.
997 RooArgSet plotClones;
998 for(std::size_t index= 0; index < plotVars.size(); index++) {
999 const RooAbsArg *var= plotVars.at(index);
1000 const RooRealVar *realVar= dynamic_cast<const RooRealVar*>(var);
1001 if(nullptr == realVar) {
1002 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot plot variable \"" << var->GetName()
1003 << "\" of type " << var->ClassName() << std::endl;
1004 return nullptr;
1005 }
1006 if(!this->dependsOn(*realVar)) {
1007 coutE(InputArguments) << ClassName() << "::" << GetName()
1008 << ":fillHistogram: WARNING: variable is not an explicit dependent: " << realVar->GetName() << std::endl;
1009 }
1010 plotClones.addClone(*realVar,true); // do not complain about duplicates
1011 }
1012
1013 // Reconnect all plotClones to each other, imported when plotting N-dim integrals with entangled parameterized ranges
1014 for(RooAbsArg * pc : plotClones) {
1015 pc->recursiveRedirectServers(plotClones,false,false,true) ;
1016 }
1017
1018 // Call checkObservables
1019 RooArgSet allDeps(plotClones) ;
1020 if (projectedVars) {
1021 allDeps.add(*projectedVars) ;
1022 }
1023 if (checkObservables(&allDeps)) {
1024 coutE(InputArguments) << "RooAbsReal::fillHistogram(" << GetName() << ") error in checkObservables, abort" << std::endl ;
1025 return hist ;
1026 }
1027
1028 // Create a standalone projection object to use for calculating bin contents
1029 RooArgSet *cloneSet = nullptr;
1030 const RooAbsReal *projected= createPlotProjection(plotClones,projectedVars,cloneSet,nullptr,condObs);
1031
1032 cxcoutD(Plotting) << "RooAbsReal::fillHistogram(" << GetName() << ") plot projection object is " << projected->GetName() << std::endl ;
1033
1034 // Prepare to loop over the histogram bins
1035 Int_t xbins(0);
1036 Int_t ybins(1);
1037 Int_t zbins(1);
1038 RooRealVar *xvar = nullptr;
1039 RooRealVar *yvar = nullptr;
1040 RooRealVar *zvar = nullptr;
1041 TAxis *xaxis = nullptr;
1042 TAxis *yaxis = nullptr;
1043 TAxis *zaxis = nullptr;
1044 switch(hdim) {
1045 case 3:
1046 zbins= hist->GetNbinsZ();
1047 zvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(2)->GetName()));
1048 zaxis= hist->GetZaxis();
1049 assert(nullptr != zvar && nullptr != zaxis);
1050 if (scaleForDensity) {
1051 scaleFactor*= (zaxis->GetXmax() - zaxis->GetXmin())/zbins;
1052 }
1053 // fall through to next case...
1054 case 2:
1055 ybins= hist->GetNbinsY();
1056 yvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(1)->GetName()));
1057 yaxis= hist->GetYaxis();
1058 assert(nullptr != yvar && nullptr != yaxis);
1059 if (scaleForDensity) {
1060 scaleFactor*= (yaxis->GetXmax() - yaxis->GetXmin())/ybins;
1061 }
1062 // fall through to next case...
1063 case 1:
1064 xbins= hist->GetNbinsX();
1065 xvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(0)->GetName()));
1066 xaxis= hist->GetXaxis();
1067 assert(nullptr != xvar && nullptr != xaxis);
1068 if (scaleForDensity) {
1069 scaleFactor*= (xaxis->GetXmax() - xaxis->GetXmin())/xbins;
1070 }
1071 break;
1072 default:
1073 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot fill histogram with "
1074 << hdim << " dimensions" << std::endl;
1075 break;
1076 }
1077
1078 // Loop over the input histogram's bins and fill each one with our projection's
1079 // value, calculated at the center.
1081 Int_t xbin(0);
1082 Int_t ybin(0);
1083 Int_t zbin(0);
1084 Int_t bins= xbins*ybins*zbins;
1085 for(Int_t bin= 0; bin < bins; bin++) {
1086 switch(hdim) {
1087 case 3:
1088 if(bin % (xbins*ybins) == 0) {
1089 zbin++;
1090 zvar->setVal(zaxis->GetBinCenter(zbin));
1091 }
1092 // fall through to next case...
1093 case 2:
1094 if(bin % xbins == 0) {
1095 ybin= (ybin%ybins) + 1;
1096 yvar->setVal(yaxis->GetBinCenter(ybin));
1097 }
1098 // fall through to next case...
1099 case 1:
1100 xbin= (xbin%xbins) + 1;
1101 xvar->setVal(xaxis->GetBinCenter(xbin));
1102 break;
1103 default:
1104 coutE(InputArguments) << "RooAbsReal::fillHistogram: Internal Error!" << std::endl;
1105 break;
1106 }
1107
1108 double result= scaleFactor*projected->getVal();
1109 if (RooAbsReal::numEvalErrors()>0) {
1110 coutW(Plotting) << "WARNING: Function evaluation error(s) at coordinates [x]=" << xvar->getVal() ;
1111 if (hdim==2) ccoutW(Plotting) << " [y]=" << yvar->getVal() ;
1112 if (hdim==3) ccoutW(Plotting) << " [z]=" << zvar->getVal() ;
1113 ccoutW(Plotting) << std::endl ;
1114 // RooAbsReal::printEvalErrors(ccoutW(Plotting),10) ;
1115 result = 0 ;
1116 }
1118
1119 hist->SetBinContent(hist->GetBin(xbin,ybin,zbin),result);
1120 if (setError) {
1121 hist->SetBinError(hist->GetBin(xbin,ybin,zbin),sqrt(result)) ;
1122 }
1123
1124 //cout << "bin " << bin << " -> (" << xbin << "," << ybin << "," << zbin << ") = " << result << std::endl;
1125 }
1127
1128 // cleanup
1129 delete cloneSet;
1130
1131 return hist;
1132}
1133
1134
1135
1136////////////////////////////////////////////////////////////////////////////////
1137/// Fill a RooDataHist with values sampled from this function at the
1138/// bin centers. If extendedMode is true, the p.d.f. values is multiplied
1139/// by the number of expected events in each bin
1140///
1141/// An optional scaling by a given scaleFactor can be performed.
1142/// Returns a pointer to the input RooDataHist, or zero
1143/// in case of an error.
1144///
1145/// If correctForBinSize is true the RooDataHist
1146/// is filled with the functions density (function value times the
1147/// bin volume) rather than function value.
1148///
1149/// If showProgress is true
1150/// a process indicator is printed on stdout in steps of one percent,
1151/// which is mostly useful for the sampling of expensive functions
1152/// such as likelihoods
1153
1154RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSet, double scaleFactor,
1155 bool correctForBinSize, bool showProgress) const
1156{
1157 // Do we have a valid histogram to use?
1158 if(nullptr == hist) {
1159 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillDataHist: no valid RooDataHist to fill" << std::endl;
1160 return nullptr;
1161 }
1162
1163 // Call checkObservables
1164 RooArgSet allDeps(*hist->get()) ;
1165 if (checkObservables(&allDeps)) {
1166 coutE(InputArguments) << "RooAbsReal::fillDataHist(" << GetName() << ") error in checkObservables, abort" << std::endl ;
1167 return hist ;
1168 }
1169
1170 // Make deep clone of self and attach to dataset observables
1171 //RooArgSet* origObs = getObservables(hist) ;
1172 RooArgSet cloneSet;
1173 RooArgSet(*this).snapshot(cloneSet, true);
1174 RooAbsReal* theClone = static_cast<RooAbsReal*>(cloneSet.find(GetName()));
1175 theClone->recursiveRedirectServers(*hist->get()) ;
1176 //const_cast<RooAbsReal*>(this)->recursiveRedirectServers(*hist->get()) ;
1177
1178 // Iterator over all bins of RooDataHist and fill weights
1179 Int_t onePct = hist->numEntries()/100 ;
1180 if (onePct==0) {
1181 onePct++ ;
1182 }
1183 for (Int_t i=0 ; i<hist->numEntries() ; i++) {
1184 if (showProgress && (i%onePct==0)) {
1185 ccoutP(Eval) << "." << std::flush ;
1186 }
1187 const RooArgSet* obs = hist->get(i) ;
1188 double binVal = theClone->getVal(normSet?normSet:obs)*scaleFactor ;
1189 if (correctForBinSize) {
1190 binVal*= hist->binVolume() ;
1191 }
1192 hist->set(i, binVal, 0.);
1193 }
1194
1195 return hist;
1196}
1197
1198
1199
1200
1201////////////////////////////////////////////////////////////////////////////////
1202/// Create and fill a ROOT histogram TH1, TH2 or TH3 with the values of this function for the variables with given names.
1203/// \param[in] varNameList List of variables to use for x, y, z axis, separated by ':'
1204/// \param[in] xbins Number of bins for first variable
1205/// \param[in] ybins Number of bins for second variable
1206/// \param[in] zbins Number of bins for third variable
1207/// \return TH1*, which is one of TH[1-3]. The histogram is owned by the caller.
1208///
1209/// For a greater degree of control use
1210/// RooAbsReal::createHistogram(const char *, const RooAbsRealLValue&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&) const
1211///
1212
1213TH1* RooAbsReal::createHistogram(RooStringView varNameList, Int_t xbins, Int_t ybins, Int_t zbins) const
1214{
1215 std::unique_ptr<RooArgSet> vars{getVariables()};
1216
1217 auto varNames = ROOT::Split(varNameList, ",:");
1218 std::vector<RooRealVar*> histVars(3, nullptr);
1219
1220 for(std::size_t iVar = 0; iVar < varNames.size(); ++iVar) {
1221 if(varNames[iVar].empty()) continue;
1222 if(iVar >= 3) {
1223 std::stringstream errMsg;
1224 errMsg << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR more than three variable names passed, but maximum number of supported variables is three";
1225 coutE(Plotting) << errMsg.str() << std::endl;
1226 throw std::invalid_argument(errMsg.str());
1227 }
1228 auto var = static_cast<RooRealVar*>(vars->find(varNames[iVar].c_str()));
1229 if(!var) {
1230 std::stringstream errMsg;
1231 errMsg << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR variable " << varNames[iVar] << " does not exist in argset: " << *vars;
1232 coutE(Plotting) << errMsg.str() << std::endl;
1233 throw std::runtime_error(errMsg.str());
1234 }
1235 histVars[iVar] = var;
1236 }
1237
1238 // Construct list of named arguments to pass to the implementation version of createHistogram()
1239
1240 RooLinkedList argList ;
1241 if (xbins>0) {
1242 argList.Add(RooFit::Binning(xbins).Clone()) ;
1243 }
1244
1245 if (histVars[1]) {
1246 argList.Add(RooFit::YVar(*histVars[1], ybins > 0 ? RooFit::Binning(ybins) : RooCmdArg::none()).Clone()) ;
1247 }
1248
1249 if (histVars[2]) {
1250 argList.Add(RooFit::ZVar(*histVars[2], zbins > 0 ? RooFit::Binning(zbins) : RooCmdArg::none()).Clone()) ;
1251 }
1252
1253 // Call implementation function
1254 TH1* result = createHistogram(GetName(), *histVars[0], argList) ;
1255
1256 // Delete temporary list of RooCmdArgs
1257 argList.Delete() ;
1258
1259 return result ;
1260}
1261
1262
1263
1264////////////////////////////////////////////////////////////////////////////////
1265/// Create and fill a ROOT histogram TH1, TH2 or TH3 with the values of this function.
1266///
1267/// \param[in] name Name of the ROOT histogram
1268/// \param[in] xvar Observable to be std::mapped on x axis of ROOT histogram
1269/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 Arguments according to list below
1270/// \return TH1 *, one of TH{1,2,3}. The caller takes ownership.
1271///
1272/// <table>
1273/// <tr><th><th> Effect on histogram creation
1274/// <tr><td> `IntrinsicBinning()` <td> Apply binning defined by function or pdf (as advertised via binBoundaries() method)
1275/// <tr><td> `Binning(const char* name)` <td> Apply binning with given name to x axis of histogram
1276/// <tr><td> `Binning(RooAbsBinning& binning)` <td> Apply specified binning to x axis of histogram
1277/// <tr><td> `Binning(int nbins, [double lo, double hi])` <td> Apply specified binning to x axis of histogram
1278/// <tr><td> `ConditionalObservables(Args_t &&... argsOrArgSet)` <td> Do not normalise PDF over following observables when projecting PDF into histogram.
1279// Arguments can either be multiple RooRealVar or a single RooArgSet containing them.
1280/// <tr><td> `Scaling(bool)` <td> Apply density-correction scaling (multiply by bin volume), default is true
1281/// <tr><td> `Extended(bool)` <td> Plot event yield instead of probability density (for extended pdfs only)
1282///
1283/// <tr><td> `YVar(const RooAbsRealLValue& var,...)` <td> Observable to be std::mapped on y axis of ROOT histogram.
1284/// The YVar() and ZVar() arguments can be supplied with optional Binning() arguments to control the binning of the Y and Z axes, e.g.
1285/// ```
1286/// createHistogram("histo",x,Binning(-1,1,20), YVar(y,Binning(-1,1,30)), ZVar(z,Binning("zbinning")))
1287/// ```
1288/// <tr><td> `ZVar(const RooAbsRealLValue& var,...)` <td> Observable to be std::mapped on z axis of ROOT histogram
1289/// </table>
1290///
1291///
1292
1294 const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
1295 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
1296{
1297
1299 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
1300 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
1301 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
1302 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
1303
1304 return createHistogram(name,xvar,l) ;
1305}
1306
1307
1308////////////////////////////////////////////////////////////////////////////////
1309/// Internal method implementing createHistogram
1310
1311TH1* RooAbsReal::createHistogram(const char *name, const RooAbsRealLValue& xvar, RooLinkedList& argList) const
1312{
1313
1314 // Define configuration for this method
1315 RooCmdConfig pc("RooAbsReal::createHistogram(" + std::string(GetName()) + ")");
1316 pc.defineInt("scaling","Scaling",0,1) ;
1317 pc.defineInt("intBinning","IntrinsicBinning",0,2) ;
1318 pc.defineInt("extended","Extended",0,2) ;
1319
1320 pc.defineSet("compSet","SelectCompSet",0);
1321 pc.defineString("compSpec","SelectCompSpec",0) ;
1322 pc.defineSet("projObs","ProjectedObservables",0,nullptr) ;
1323 pc.defineObject("yvar","YVar",0,nullptr) ;
1324 pc.defineObject("zvar","ZVar",0,nullptr) ;
1325 pc.defineMutex("SelectCompSet","SelectCompSpec") ;
1326 pc.defineMutex("IntrinsicBinning","Binning") ;
1327 pc.defineMutex("IntrinsicBinning","BinningName") ;
1328 pc.defineMutex("IntrinsicBinning","BinningSpec") ;
1329 pc.allowUndefined() ;
1330
1331 // Process & check varargs
1332 pc.process(argList) ;
1333 if (!pc.ok(true)) {
1334 return nullptr ;
1335 }
1336
1337 RooArgList vars(xvar) ;
1338 RooAbsArg* yvar = static_cast<RooAbsArg*>(pc.getObject("yvar")) ;
1339 if (yvar) {
1340 vars.add(*yvar) ;
1341 }
1342 RooAbsArg* zvar = static_cast<RooAbsArg*>(pc.getObject("zvar")) ;
1343 if (zvar) {
1344 vars.add(*zvar) ;
1345 }
1346
1347 auto projObs = pc.getSet("projObs");
1348 RooArgSet* intObs = nullptr ;
1349
1350 bool doScaling = pc.getInt("scaling") ;
1351 Int_t doIntBinning = pc.getInt("intBinning") ;
1352 Int_t doExtended = pc.getInt("extended") ;
1353
1354 // If doExtended is two, selection is automatic, set to 1 of pdf is extended, to zero otherwise
1355 const RooAbsPdf* pdfSelf = dynamic_cast<const RooAbsPdf*>(this) ;
1356 if (!pdfSelf && doExtended == 1) {
1357 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-pdf object, ignored" << std::endl ;
1358 doExtended=0 ;
1359 }
1360 if (pdfSelf && doExtended==1 && pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended) {
1361 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-extendable pdf, ignored" << std::endl ;
1362 doExtended=0 ;
1363 }
1364 if (pdfSelf && doExtended==2) {
1365 doExtended = pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended ? 0 : 1 ;
1366 } else if(!pdfSelf) {
1367 doExtended = 0;
1368 }
1369
1370 const char* compSpec = pc.getString("compSpec") ;
1371 const RooArgSet* compSet = pc.getSet("compSet");
1372 bool haveCompSel = ( (compSpec && strlen(compSpec)>0) || compSet) ;
1373
1374 std::unique_ptr<RooBinning> intBinning;
1375 if (doIntBinning>0) {
1376 // Given RooAbsPdf* pdf and RooRealVar* obs
1377 std::unique_ptr<std::list<double>> bl{binBoundaries(const_cast<RooAbsRealLValue&>(xvar),xvar.getMin(),xvar.getMax())};
1378 if (!bl) {
1379 // Only emit warning when intrinsic binning is explicitly requested
1380 if (doIntBinning==1) {
1381 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
1382 << ") WARNING, intrinsic model binning requested for histogram, but model does not define bin boundaries, reverting to default binning"<< std::endl ;
1383 }
1384 } else {
1385 if (doIntBinning==2) {
1386 coutI(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
1387 << ") INFO: Model has intrinsic binning definition, selecting that binning for the histogram"<< std::endl ;
1388 }
1389 std::vector<double> edges(bl->size());
1390 int i=0 ;
1391 for (auto const& elem : *bl) { edges[i++] = elem ; }
1392 intBinning = std::make_unique<RooBinning>(bl->size()-1,edges.data()) ;
1393 }
1394 }
1395
1396 RooLinkedList argListCreate(argList) ;
1397 RooCmdConfig::stripCmdList(argListCreate,"Scaling,ProjectedObservables,IntrinsicBinning,SelectCompSet,SelectCompSpec,Extended") ;
1398
1399 TH1* histo(nullptr) ;
1400 if (intBinning) {
1401 RooCmdArg tmp = RooFit::Binning(*intBinning) ;
1402 argListCreate.Add(&tmp) ;
1403 histo = xvar.createHistogram(name,argListCreate) ;
1404 } else {
1405 histo = xvar.createHistogram(name,argListCreate) ;
1406 }
1407
1408 // Do component selection here
1409 if (haveCompSel) {
1410
1411 // Get complete set of tree branch nodes
1412 RooArgSet branchNodeSet ;
1413 branchNodeServerList(&branchNodeSet) ;
1414
1415 // Discard any non-RooAbsReal nodes
1416 for(RooAbsArg * arg : branchNodeSet) {
1417 if (!dynamic_cast<RooAbsReal*>(arg)) {
1418 branchNodeSet.remove(*arg) ;
1419 }
1420 }
1421
1422 std::unique_ptr<RooArgSet> dirSelNodes;
1423 if (compSet) {
1424 dirSelNodes.reset(static_cast<RooArgSet*>(branchNodeSet.selectCommon(*compSet)));
1425 } else {
1426 dirSelNodes.reset(static_cast<RooArgSet*>(branchNodeSet.selectByName(compSpec)));
1427 }
1428 if (!dirSelNodes->empty()) {
1429 coutI(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") directly selected PDF components: " << *dirSelNodes << std::endl ;
1430
1431 // Do indirect selection and activate both
1432 plotOnCompSelect(dirSelNodes.get()) ;
1433 } else {
1434 if (compSet) {
1435 coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection set " << *compSet << " does not match any components of p.d.f." << std::endl ;
1436 } else {
1437 coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection expression '" << compSpec << "' does not select any components of p.d.f." << std::endl ;
1438 }
1439 return nullptr ;
1440 }
1441 }
1442
1443 double scaleFactor(1.0) ;
1444 if (doExtended) {
1445 scaleFactor = pdfSelf->expectedEvents(vars) ;
1446 doScaling=false ;
1447 }
1448
1449 fillHistogram(histo,vars,scaleFactor,intObs,doScaling,projObs,false) ;
1450
1451 // Deactivate component selection
1452 if (haveCompSel) {
1453 plotOnCompSelect(nullptr) ;
1454 }
1455
1456
1457 return histo ;
1458}
1459
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Helper function for plotting of composite p.d.fs. Given
1463/// a set of selected components that should be plotted,
1464/// find all nodes that (in)directly depend on these selected
1465/// nodes. Mark all directly and indirectly selected nodes
1466/// as 'selected' using the selectComp() method
1467
1469{
1470 // Get complete set of tree branch nodes
1471 RooArgSet branchNodeSet;
1472 branchNodeServerList(&branchNodeSet);
1473
1474 // Discard any non-PDF nodes
1475 // Iterate by number because collection is being modified! Iterators may invalidate ...
1476 for (unsigned int i = 0; i < branchNodeSet.size(); ++i) {
1477 const auto arg = branchNodeSet[i];
1478 if (!dynamic_cast<RooAbsReal*>(arg)) {
1479 branchNodeSet.remove(*arg) ;
1480 }
1481 }
1482
1483 // If no set is specified, restored all selection bits to true
1484 if (!selNodes) {
1485 // Reset PDF selection bits to true
1486 for (const auto arg : branchNodeSet) {
1487 static_cast<RooAbsReal*>(arg)->selectComp(true);
1488 }
1489 return ;
1490 }
1491
1492
1493 // Add all nodes below selected nodes that are value servers
1494 RooArgSet tmp;
1495 for (const auto arg : branchNodeSet) {
1496 for (const auto selNode : *selNodes) {
1497 if (selNode->dependsOn(*arg, nullptr, /*valueOnly=*/true)) {
1498 tmp.add(*arg,true);
1499 }
1500 }
1501 }
1502
1503 // Add all nodes that depend on selected nodes by value
1504 for (const auto arg : branchNodeSet) {
1505 if (arg->dependsOn(*selNodes, nullptr, /*valueOnly=*/true)) {
1506 tmp.add(*arg,true);
1507 }
1508 }
1509
1510 tmp.remove(*selNodes, true);
1511 tmp.remove(*this);
1512 selNodes->add(tmp);
1513 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") indirectly selected PDF components: " << tmp << std::endl ;
1514
1515 // Set PDF selection bits according to selNodes
1516 for (const auto arg : branchNodeSet) {
1517 bool select = selNodes->find(arg->GetName()) != nullptr;
1518 static_cast<RooAbsReal*>(arg)->selectComp(select);
1519 }
1520}
1521
1522
1523
1524////////////////////////////////////////////////////////////////////////////////
1525/// Plot (project) PDF on specified frame. If a PDF is plotted in an empty frame, it
1526/// will show a unit normalized curve in the frame variable, taken at the present value
1527/// of other observables defined for this PDF.
1528///
1529/// \param[in] frame pointer to RooPlot
1530/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10 Ordered arguments
1531///
1532/// If a PDF is plotted in a frame in which a dataset has already been plotted, it will
1533/// show a projected curve integrated over all variables that were present in the shown
1534/// dataset except for the one on the x-axis. The normalization of the curve will also
1535/// be adjusted to the event count of the plotted dataset. An informational message
1536/// will be printed for each projection step that is performed.
1537///
1538/// This function takes the following named arguments
1539/// <table>
1540/// <tr><th><th> Projection control
1541/// <tr><td> `Slice(const RooArgSet& set)` <td> Override default projection behaviour by omitting observables listed
1542/// in set from the projection, i.e. by not integrating over these.
1543/// Slicing is usually only sensible in discrete observables, by e.g. creating a slice
1544/// of the PDF at the current value of the category observable.
1545///
1546/// <tr><td> `Slice(RooCategory& cat, const char* label)` <td> Override default projection behaviour by omitting the specified category
1547/// observable from the projection, i.e., by not integrating over all states of this category.
1548/// The slice is positioned at the given label value. To pass multiple Slice() commands, please use the
1549/// Slice(std::map<RooCategory*, std::string> const&) argument explained below.
1550///
1551/// <tr><td> `Slice(std::map<RooCategory*, std::string> const&)` <td> Omits multiple categories from the projection, as explianed above.
1552/// Can be used with initializer lists for convenience, e.g.
1553/// ```{.cpp}
1554/// pdf.plotOn(frame, Slice({{&tagCategory, "2tag"}, {&jetCategory, "3jet"}});
1555/// ```
1556///
1557/// <tr><td> `Project(const RooArgSet& set)` <td> Override default projection behaviour by projecting over observables
1558/// given in the set, ignoring the default projection behavior. Advanced use only.
1559///
1560/// <tr><td> `ProjWData(const RooAbsData& d)` <td> Override default projection _technique_ (integration). For observables present in given dataset
1561/// projection of PDF is achieved by constructing an average over all observable values in given set.
1562/// Consult RooFit plotting tutorial for further explanation of meaning & use of this technique
1563///
1564/// <tr><td> `ProjWData(const RooArgSet& s, const RooAbsData& d)` <td> As above but only consider subset 's' of observables in dataset 'd' for projection through data averaging
1565///
1566/// <tr><td> `ProjectionRange(const char* rn)` <td> Override default range of projection integrals to a different range specified by given range name.
1567/// This technique allows you to project a finite width slice in a real-valued observable
1568///
1569/// <tr><td> `NumCPU(Int_t ncpu)` <td> Number of CPUs to use simultaneously to calculate data-weighted projections (only in combination with ProjWData)
1570///
1571///
1572/// <tr><th><th> Misc content control
1573/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per curve. A negative
1574/// value suppress output completely, a zero value will only print the error count per p.d.f component,
1575/// a positive value is will print details of each error up to numErr messages per p.d.f component.
1576///
1577/// <tr><td> `EvalErrorValue(double value)` <td> Set curve points at which (pdf) evaluation errors occur to specified value. By default the
1578/// function value is plotted.
1579///
1580/// <tr><td> `Normalization(double scale, ScaleType code)` <td> Adjust normalization by given scale factor. Interpretation of number depends on code:
1581/// - Relative: relative adjustment factor for a normalized function,
1582/// - NumEvent: scale to match given number of events.
1583/// - Raw: relative adjustment factor for an un-normalized function.
1584///
1585/// <tr><td> `Name(const chat* name)` <td> Give curve specified name in frame. Useful if curve is to be referenced later
1586///
1587/// <tr><td> `Asymmetry(const RooCategory& c)` <td> Show the asymmetry of the PDF in given two-state category [F(+)-F(-)] / [F(+)+F(-)] rather than
1588/// the PDF projection. Category must have two states with indices -1 and +1 or three states with
1589/// indices -1,0 and +1.
1590///
1591/// <tr><td> `ShiftToZero(bool flag)` <td> Shift entire curve such that lowest visible point is at exactly zero. Mostly useful when plotting \f$ -\log(L) \f$ or \f$ \chi^2 \f$ distributions
1592///
1593/// <tr><td> `AddTo(const char* name, double_t wgtSelf, double_t wgtOther)` <td> Add constructed projection to already existing curve with given name and relative weight factors
1594/// <tr><td> `Components(const char* names)` <td> When plotting sums of PDFs, plot only the named components (*e.g.* only
1595/// the signal of a signal+background model).
1596/// <tr><td> `Components(const RooArgSet& compSet)` <td> As above, but pass a RooArgSet of the components themselves.
1597///
1598/// <tr><th><th> Plotting control
1599/// <tr><td> `DrawOption(const char* opt)` <td> Select ROOT draw option for resulting TGraph object. Currently supported options are "F" (fill), "L" (line), and "P" (points).
1600/// \note Option "P" will cause RooFit to plot (and treat) this pdf as if it were data! This is intended for plotting "corrected data"-type pdfs such as "data-minus-background" or unfolded datasets.
1601///
1602/// <tr><td> `LineStyle(Int_t style)` <td> Select line style by ROOT line style code, default is solid
1603///
1604/// <tr><td> `LineColor(Int_t color)` <td> Select line color by ROOT color code, default is blue
1605///
1606/// <tr><td> `LineWidth(Int_t width)` <td> Select line with in pixels, default is 3
1607///
1608/// <tr><td> `MarkerStyle(Int_t style)` <td> Select the ROOT marker style, default is 21
1609///
1610/// <tr><td> `MarkerColor(Int_t color)` <td> Select the ROOT marker color, default is black
1611///
1612/// <tr><td> `MarkerSize(double size)` <td> Select the ROOT marker size
1613///
1614/// <tr><td> `FillStyle(Int_t style)` <td> Select fill style, default is not filled. If a filled style is selected, also use VLines()
1615/// to add vertical downward lines at end of curve to ensure proper closure. Add `DrawOption("F")` for filled drawing.
1616/// <tr><td> `FillColor(Int_t color)` <td> Select fill color by ROOT color code
1617///
1618/// <tr><td> `Range(const char* name)` <td> Only draw curve in range defined by given name
1619///
1620/// <tr><td> `Range(double lo, double hi)` <td> Only draw curve in specified range
1621///
1622/// <tr><td> `VLines()` <td> Add vertical lines to y=0 at end points of curve
1623///
1624/// <tr><td> `Precision(double eps)` <td> Control precision of drawn curve w.r.t to scale of plot, default is 1e-3. Higher precision
1625/// will result in more and more densely spaced curve points
1626///
1627/// <tr><td> `Invisible(bool flag)` <td> Add curve to frame, but do not display. Useful in combination AddTo()
1628///
1629/// <tr><td> `VisualizeError(const RooFitResult& fitres, double Z=1, bool linearMethod=true)`
1630/// <td> Visualize the uncertainty on the parameters, as given in fitres, at 'Z' sigma'. The linear method is fast but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made. 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
1631///
1632/// <tr><td> `VisualizeError(const RooFitResult& fitres, const RooArgSet& param, double Z=1, bool linearMethod=true)`
1633/// <td> Visualize the uncertainty on the subset of parameters 'param', as given in fitres, at 'Z' sigma'
1634/// </table>
1635///
1636/// Details on error band visualization
1637/// -----------------------------------
1638/// *VisualizeError() uses plotOnWithErrorBand(). Documentation of the latter:*
1639/// \see plotOnWithErrorBand()
1640
1641RooPlot* RooAbsReal::plotOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2,
1642 const RooCmdArg& arg3, const RooCmdArg& arg4,
1643 const RooCmdArg& arg5, const RooCmdArg& arg6,
1644 const RooCmdArg& arg7, const RooCmdArg& arg8,
1645 const RooCmdArg& arg9, const RooCmdArg& arg10) const
1646{
1648 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
1649 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
1650 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
1651 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
1652 l.Add((TObject*)&arg9) ; l.Add((TObject*)&arg10) ;
1653 return plotOn(frame,l) ;
1654}
1655
1656
1657
1658////////////////////////////////////////////////////////////////////////////////
1659/// Internal back-end function of plotOn() with named arguments
1660
1662{
1663 // Special handling here if argList contains RangeWithName argument with multiple
1664 // range names -- Need to translate this call into multiple calls
1665
1666 RooCmdArg* rcmd = static_cast<RooCmdArg*>(argList.FindObject("RangeWithName")) ;
1667 if (rcmd && TString(rcmd->getString(0)).Contains(",")) {
1668
1669 // List joint ranges as choice of normalization for all later processing
1670 RooCmdArg rnorm = RooFit::NormRange(rcmd->getString(0)) ;
1671 argList.Add(&rnorm) ;
1672
1673 std::vector<std::string> rlist;
1674
1675 // Separate named ranges using strtok
1676 for (const std::string& rangeNameToken : ROOT::Split(rcmd->getString(0), ",")) {
1677 rlist.emplace_back(rangeNameToken);
1678 }
1679
1680 for (const auto& rangeString : rlist) {
1681 // Process each range with a separate command with a single range to be plotted
1682 rcmd->setString(0, rangeString.c_str());
1683 RooAbsReal::plotOn(frame,argList);
1684 }
1685 return frame ;
1686
1687 }
1688
1689 // Define configuration for this method
1690 RooCmdConfig pc("RooAbsReal::plotOn(" + std::string(GetName()) + ")");
1691 pc.defineString("drawOption","DrawOption",0,"L") ;
1692 pc.defineString("projectionRangeName","ProjectionRange",0,"",true) ;
1693 pc.defineString("curveNameSuffix","CurveNameSuffix",0,"") ;
1694 pc.defineString("sliceCatState","SliceCat",0,"",true) ;
1695 pc.defineDouble("scaleFactor","Normalization",0,1.0) ;
1696 pc.defineInt("scaleType","Normalization",0,Relative) ;
1697 pc.defineSet("sliceSet","SliceVars",0) ;
1698 pc.defineObject("sliceCatList","SliceCat",0,nullptr,true) ;
1699 // This dummy is needed for plotOn to recognize the "SliceCatMany" command.
1700 // It is not used directly, but the "SliceCat" commands are nested in it.
1701 // Removing this dummy definition results in "ERROR: unrecognized command: SliceCatMany".
1702 pc.defineObject("dummy1","SliceCatMany",0) ;
1703 pc.defineSet("projSet","Project",0) ;
1704 pc.defineObject("asymCat","Asymmetry",0) ;
1705 pc.defineDouble("precision","Precision",0,1e-3) ;
1706 pc.defineDouble("evalErrorVal","EvalErrorValue",0,0) ;
1707 pc.defineInt("doEvalError","EvalErrorValue",0,0) ;
1708 pc.defineInt("shiftToZero","ShiftToZero",0,0) ;
1709 pc.defineSet("projDataSet","ProjData",0) ;
1710 pc.defineObject("projData","ProjData",1) ;
1711 pc.defineObject("errorFR","VisualizeError",0) ;
1712 pc.defineDouble("errorZ","VisualizeError",0,1.) ;
1713 pc.defineSet("errorPars","VisualizeError",0) ;
1714 pc.defineInt("linearMethod","VisualizeError",0,0) ;
1715 pc.defineInt("binProjData","ProjData",0,0) ;
1716 pc.defineDouble("rangeLo","Range",0,-999.) ;
1717 pc.defineDouble("rangeHi","Range",1,-999.) ;
1718 pc.defineInt("numee","PrintEvalErrors",0,10) ;
1719 pc.defineInt("rangeAdjustNorm","Range",0,0) ;
1720 pc.defineInt("rangeWNAdjustNorm","RangeWithName",0,0) ;
1721 pc.defineInt("VLines","VLines",0,2) ; // 2==ExtendedWings
1722 pc.defineString("rangeName","RangeWithName",0,"") ;
1723 pc.defineString("normRangeName","NormRange",0,"") ;
1724 pc.defineInt("markerColor","MarkerColor",0,-999) ;
1725 pc.defineInt("markerStyle","MarkerStyle",0,-999) ;
1726 pc.defineDouble("markerSize","MarkerSize",0,-999) ;
1727 pc.defineInt("lineColor","LineColor",0,-999) ;
1728 pc.defineInt("lineStyle","LineStyle",0,-999) ;
1729 pc.defineInt("lineWidth","LineWidth",0,-999) ;
1730 pc.defineInt("fillColor","FillColor",0,-999) ;
1731 pc.defineInt("fillStyle","FillStyle",0,-999) ;
1732 pc.defineString("curveName","Name",0,"") ;
1733 pc.defineInt("curveInvisible","Invisible",0,0) ;
1734 pc.defineInt("showProg","ShowProgress",0,0) ;
1735 pc.defineInt("numCPU","NumCPU",0,1) ;
1736 pc.defineInt("interleave","NumCPU",1,0) ;
1737 pc.defineString("addToCurveName","AddTo",0,"") ;
1738 pc.defineDouble("addToWgtSelf","AddTo",0,1.) ;
1739 pc.defineDouble("addToWgtOther","AddTo",1,1.) ;
1740 pc.defineInt("moveToBack","MoveToBack",0,0) ;
1741 pc.defineMutex("SliceVars","Project") ;
1742 pc.defineMutex("AddTo","Asymmetry") ;
1743 pc.defineMutex("Range","RangeWithName") ;
1744 pc.defineMutex("VisualizeError","VisualizeErrorData") ;
1745
1746 // Process & check varargs
1747 pc.process(argList) ;
1748 if (!pc.ok(true)) {
1749 return frame ;
1750 }
1751
1752 PlotOpt o ;
1753 TString drawOpt(pc.getString("drawOption"));
1754
1755 RooFitResult* errFR = static_cast<RooFitResult*>(pc.getObject("errorFR")) ;
1756 double errZ = pc.getDouble("errorZ") ;
1757 RooArgSet* errPars = pc.getSet("errorPars") ;
1758 bool linMethod = pc.getInt("linearMethod") ;
1759 if (!drawOpt.Contains("P") && errFR) {
1760 return plotOnWithErrorBand(frame,*errFR,errZ,errPars,argList,linMethod) ;
1761 } else {
1762 o.errorFR = errFR;
1763 }
1764
1765 // Extract values from named arguments
1766 o.numee = pc.getInt("numee") ;
1767 o.drawOptions = drawOpt.Data();
1768 o.curveNameSuffix = pc.getString("curveNameSuffix") ;
1769 o.scaleFactor = pc.getDouble("scaleFactor") ;
1770 o.stype = (ScaleType) pc.getInt("scaleType") ;
1771 o.projData = static_cast<const RooAbsData*>(pc.getObject("projData")) ;
1772 o.binProjData = pc.getInt("binProjData") ;
1773 o.projDataSet = pc.getSet("projDataSet");
1774 o.numCPU = pc.getInt("numCPU") ;
1775 o.interleave = (RooFit::MPSplit) pc.getInt("interleave") ;
1776 o.eeval = pc.getDouble("evalErrorVal") ;
1777 o.doeeval = pc.getInt("doEvalError") ;
1778
1779 const RooArgSet* sliceSetTmp = pc.getSet("sliceSet");
1780 std::unique_ptr<RooArgSet> sliceSet{sliceSetTmp ? static_cast<RooArgSet*>(sliceSetTmp->Clone()) : nullptr};
1781 const RooArgSet* projSet = pc.getSet("projSet") ;
1782 const RooAbsCategoryLValue* asymCat = static_cast<const RooAbsCategoryLValue*>(pc.getObject("asymCat")) ;
1783
1784
1785 // Look for category slice arguments and add them to the master slice list if found
1786 const char* sliceCatState = pc.getString("sliceCatState",nullptr,true) ;
1787 const RooLinkedList& sliceCatList = pc.getObjectList("sliceCatList") ;
1788 if (sliceCatState) {
1789
1790 // Make the master slice set if it doesnt exist
1791 if (!sliceSet) {
1792 sliceSet = std::make_unique<RooArgSet>();
1793 }
1794
1795 // Loop over all categories provided by (multiple) Slice() arguments
1796 auto catTokens = ROOT::Split(sliceCatState, ",");
1797 auto iter = sliceCatList.begin();
1798 for (unsigned int i=0; i < catTokens.size(); ++i) {
1799 if (auto scat = static_cast<RooCategory*>(*iter)) {
1800 // Set the slice position to the value indicate by slabel
1801 scat->setLabel(catTokens[i]) ;
1802 // Add the slice category to the master slice set
1803 sliceSet->add(*scat,false) ;
1804 }
1805 ++iter;
1806 }
1807 }
1808
1809 o.precision = pc.getDouble("precision") ;
1810 o.shiftToZero = (pc.getInt("shiftToZero")!=0) ;
1811 Int_t vlines = pc.getInt("VLines");
1812 if (pc.hasProcessed("Range")) {
1813 o.rangeLo = pc.getDouble("rangeLo") ;
1814 o.rangeHi = pc.getDouble("rangeHi") ;
1815 o.postRangeFracScale = pc.getInt("rangeAdjustNorm") ;
1816 if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
1817 } else if (pc.hasProcessed("RangeWithName")) {
1818 o.normRangeName = pc.getString("rangeName",nullptr,true) ;
1819 o.rangeLo = frame->getPlotVar()->getMin(pc.getString("rangeName",nullptr,true)) ;
1820 o.rangeHi = frame->getPlotVar()->getMax(pc.getString("rangeName",nullptr,true)) ;
1821 o.postRangeFracScale = pc.getInt("rangeWNAdjustNorm") ;
1822 if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
1823 }
1824
1825
1826 // If separate normalization range was specified this overrides previous settings
1827 if (pc.hasProcessed("NormRange")) {
1828 o.normRangeName = pc.getString("normRangeName") ;
1829 o.postRangeFracScale = true ;
1830 }
1831
1832 o.wmode = (vlines==2)?RooCurve::Extended:(vlines==1?RooCurve::Straight:RooCurve::NoWings) ;
1833 o.projectionRangeName = pc.getString("projectionRangeName",nullptr,true) ;
1834 o.curveName = pc.getString("curveName",nullptr,true) ;
1835 o.curveInvisible = pc.getInt("curveInvisible") ;
1836 o.progress = pc.getInt("showProg") ;
1837 o.addToCurveName = pc.getString("addToCurveName",nullptr,true) ;
1838 o.addToWgtSelf = pc.getDouble("addToWgtSelf") ;
1839 o.addToWgtOther = pc.getDouble("addToWgtOther") ;
1840
1842 coutE(InputArguments) << "RooAbsReal::plotOn(" << GetName() << ") cannot find existing curve " << o.addToCurveName << " to add to in RooPlot" << std::endl ;
1843 return frame ;
1844 }
1845
1846 RooArgSet projectedVars ;
1847 if (sliceSet) {
1848 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have slice " << *sliceSet << std::endl ;
1849
1850 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1851
1852 // Take out the sliced variables
1853 for (const auto sliceArg : *sliceSet) {
1854 RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
1855 if (arg) {
1856 projectedVars.remove(*arg) ;
1857 } else {
1858 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") slice variable "
1859 << sliceArg->GetName() << " was not projected anyway" << std::endl ;
1860 }
1861 }
1862 } else if (projSet) {
1863 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have projSet " << *projSet << std::endl ;
1864 makeProjectionSet(frame->getPlotVar(),projSet,projectedVars,false) ;
1865 } else {
1866 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have neither sliceSet nor projSet " << std::endl ;
1867 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1868 }
1869 o.projSet = &projectedVars ;
1870
1871 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: projectedVars = " << projectedVars << std::endl ;
1872
1873
1874 RooPlot* ret ;
1875 if (!asymCat) {
1876 // Forward to actual calculation
1877 ret = RooAbsReal::plotOn(frame,o) ;
1878 } else {
1879 // Forward to actual calculation
1880 ret = RooAbsReal::plotAsymOn(frame,*asymCat,o) ;
1881 }
1882
1883 // Optionally adjust line/fill attributes
1884 Int_t lineColor = pc.getInt("lineColor") ;
1885 Int_t lineStyle = pc.getInt("lineStyle") ;
1886 Int_t lineWidth = pc.getInt("lineWidth") ;
1887 Int_t markerColor = pc.getInt("markerColor") ;
1888 Int_t markerStyle = pc.getInt("markerStyle") ;
1889 Size_t markerSize = pc.getDouble("markerSize") ;
1890 Int_t fillColor = pc.getInt("fillColor") ;
1891 Int_t fillStyle = pc.getInt("fillStyle") ;
1892 if (lineColor!=-999) ret->getAttLine()->SetLineColor(lineColor) ;
1893 if (lineStyle!=-999) ret->getAttLine()->SetLineStyle(lineStyle) ;
1894 if (lineWidth!=-999) ret->getAttLine()->SetLineWidth(lineWidth) ;
1895 if (fillColor!=-999) ret->getAttFill()->SetFillColor(fillColor) ;
1896 if (fillStyle!=-999) ret->getAttFill()->SetFillStyle(fillStyle) ;
1897 if (markerColor!=-999) ret->getAttMarker()->SetMarkerColor(markerColor) ;
1898 if (markerStyle!=-999) ret->getAttMarker()->SetMarkerStyle(markerStyle) ;
1899 if (markerSize!=-999) ret->getAttMarker()->SetMarkerSize(markerSize) ;
1900
1901 if ((fillColor != -999 || fillStyle != -999) && !drawOpt.Contains("F")) {
1902 coutW(Plotting) << "Fill color or style was set for plotting \"" << GetName()
1903 << "\", but these only have an effect when 'DrawOption(\"F\")' for fill is used at the same time." << std::endl;
1904 }
1905
1906 // Move last inserted object to back to drawing stack if requested
1907 if (pc.getInt("moveToBack") && frame->numItems()>1) {
1908 frame->drawBefore(frame->getObject(0)->GetName(), frame->getCurve()->GetName());
1909 }
1910
1911 return ret ;
1912}
1913
1914
1915
1916/// Plotting engine function for internal use
1917///
1918/// Plot ourselves on given frame. If frame contains a histogram, all dimensions of the plotted
1919/// function that occur in the previously plotted dataset are projected via partial integration,
1920/// otherwise no projections are performed. Optionally, certain projections can be performed
1921/// by summing over the values present in a provided dataset ('projData'), to correctly
1922/// project out data dependents that are not properly described by the PDF (e.g. per-event errors).
1923///
1924/// The functions value can be multiplied with an optional scale factor. The interpretation
1925/// of the scale factor is unique for generic real functions, for PDFs there are various interpretations
1926/// possible, which can be selection with 'stype' (see RooAbsPdf::plotOn() for details).
1927///
1928/// The default projection behaviour can be overridden by supplying an optional set of dependents
1929/// to project. For most cases, plotSliceOn() and plotProjOn() provide a more intuitive interface
1930/// to modify the default projection behaviour.
1931//_____________________________________________________________________________
1932// coverity[PASS_BY_VALUE]
1934{
1935
1936
1937 // Sanity checks
1938 if (plotSanityChecks(frame)) return frame ;
1939
1940 // ProjDataVars is either all projData observables, or the user indicated subset of it
1941 RooArgSet projDataVars ;
1942 if (o.projData) {
1943 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjData with observables = " << *o.projData->get() << std::endl ;
1944 if (o.projDataSet) {
1945 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(o.projData->get()->selectCommon(*o.projDataSet))};
1946 projDataVars.add(*tmp) ;
1947 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjDataSet = " << *o.projDataSet << " will only use this subset of projData" << std::endl ;
1948 } else {
1949 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") using full ProjData" << std::endl ;
1950 projDataVars.add(*o.projData->get()) ;
1951 }
1952 }
1953
1954 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") ProjDataVars = " << projDataVars << std::endl ;
1955
1956 // Make list of variables to be projected
1957 RooArgSet projectedVars ;
1958 RooArgSet sliceSet ;
1959 if (o.projSet) {
1960 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have input projSet = " << *o.projSet << std::endl ;
1961 makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,false) ;
1962 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") calculated projectedVars = " << *o.projSet << std::endl ;
1963
1964 // Print list of non-projected variables
1965 if (frame->getNormVars()) {
1966 RooArgSet sliceSetTmp;
1967 getObservables(frame->getNormVars(), sliceSetTmp) ;
1968
1969 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") frame->getNormVars() that are also observables = " << sliceSetTmp << std::endl ;
1970
1971 sliceSetTmp.remove(projectedVars,true,true) ;
1972 sliceSetTmp.remove(*frame->getPlotVar(),true,true) ;
1973
1974 if (o.projData) {
1975 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(projDataVars.selectCommon(*o.projSet))};
1976 sliceSetTmp.remove(*tmp,true,true) ;
1977 }
1978
1979 if (!sliceSetTmp.empty()) {
1980 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on "
1981 << frame->getPlotVar()->GetName() << " represents a slice in " << sliceSetTmp << std::endl ;
1982 }
1983 sliceSet.add(sliceSetTmp) ;
1984 }
1985 } else {
1986 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1987 }
1988
1989 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") projectedVars = " << projectedVars << " sliceSet = " << sliceSet << std::endl ;
1990
1991
1992 RooArgSet* projDataNeededVars = nullptr ;
1993 // Take out data-projected dependents from projectedVars
1994 if (o.projData) {
1995 projDataNeededVars = static_cast<RooArgSet*>(projectedVars.selectCommon(projDataVars)) ;
1996 projectedVars.remove(projDataVars,true,true) ;
1997 }
1998
1999 // Get the plot variable and remember its original value
2000 auto* plotVar = static_cast<RooRealVar*>(frame->getPlotVar());
2001 double oldPlotVarVal = plotVar->getVal();
2002
2003 // Inform user about projections
2004 if (!projectedVars.empty()) {
2005 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2006 << " integrates over variables " << projectedVars
2007 << (o.projectionRangeName?Form(" in range %s",o.projectionRangeName):"") << std::endl;
2008 }
2009 if (projDataNeededVars && !projDataNeededVars->empty()) {
2010 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2011 << " averages using data variables " << *projDataNeededVars << std::endl ;
2012 }
2013
2014 // Create projection integral
2015 RooArgSet* projectionCompList = nullptr ;
2016
2017 RooArgSet deps;
2018 getObservables(frame->getNormVars(), deps) ;
2019 deps.remove(projectedVars,true,true) ;
2020 if (projDataNeededVars) {
2021 deps.remove(*projDataNeededVars,true,true) ;
2022 }
2023 deps.remove(*plotVar,true,true) ;
2024 deps.add(*plotVar) ;
2025
2026 // Now that we have the final set of dependents, call checkObservables()
2027
2028 // WVE take out conditional observables
2029 if (checkObservables(&deps)) {
2030 coutE(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") error in checkObservables, abort" << std::endl ;
2031 if (projDataNeededVars) delete projDataNeededVars ;
2032 return frame ;
2033 }
2034
2035 RooAbsReal *projection = const_cast<RooAbsReal*>(createPlotProjection(deps, &projectedVars, projectionCompList, o.projectionRangeName));
2036 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot projection object is " << projection->GetName() << std::endl ;
2037 if (dologD(Plotting)) {
2038 projection->printStream(ccoutD(Plotting),0,kVerbose) ;
2039 }
2040
2041 // Always fix RooAddPdf normalizations
2042 RooArgSet fullNormSet(deps) ;
2043 fullNormSet.add(projectedVars) ;
2044 if (projDataNeededVars && !projDataNeededVars->empty()) {
2045 fullNormSet.add(*projDataNeededVars) ;
2046 }
2047
2048 std::unique_ptr<RooArgSet> projectionComponents(projection->getComponents());
2049 for(auto * pdf : dynamic_range_cast<RooAbsPdf*>(*projectionComponents)) {
2050 if (pdf) {
2051 pdf->selectNormalization(&fullNormSet) ;
2052 }
2053 }
2054
2055 // Apply data projection, if requested
2056 if (o.projData && projDataNeededVars && !projDataNeededVars->empty()) {
2057
2058 // If data set contains more rows than needed, make reduced copy first
2059 RooAbsData* projDataSel = const_cast<RooAbsData*>(o.projData);
2060 std::unique_ptr<RooAbsData> projDataSelOwned;
2061
2062 if (projDataNeededVars->size() < o.projData->get()->size()) {
2063
2064 // Determine if there are any slice variables in the projection set
2065 std::unique_ptr<RooArgSet> sliceDataSet{static_cast<RooArgSet*>(sliceSet.selectCommon(*o.projData->get()))};
2066 TString cutString ;
2067 if (!sliceDataSet->empty()) {
2068 bool first(true) ;
2069 for(RooAbsArg * sliceVar : *sliceDataSet) {
2070 if (!first) {
2071 cutString.Append("&&") ;
2072 } else {
2073 first=false ;
2074 }
2075
2076 RooAbsRealLValue* real ;
2078 if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
2079 cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
2080 } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
2081 cutString.Append(Form("%s==%d",cat->GetName(),cat->getCurrentIndex())) ;
2082 }
2083 }
2084 }
2085
2086 if (!cutString.IsNull()) {
2087 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") reducing given projection dataset to entries with " << cutString << std::endl ;
2088 }
2089 projDataSelOwned = std::unique_ptr<RooAbsData>{const_cast<RooAbsData*>(o.projData)->reduce(*projDataNeededVars, cutString.IsNull() ? nullptr : cutString)};
2090 projDataSel = projDataSelOwned.get();
2091 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName()
2092 << ") only the following components of the projection data will be used: " << *projDataNeededVars << std::endl ;
2093 }
2094
2095 // Request binning of unbinned projection dataset that consists exclusively of category observables
2096 if (!o.binProjData && dynamic_cast<RooDataSet*>(projDataSel)!=nullptr) {
2097
2098 // Determine if dataset contains only categories
2099 bool allCat(true) ;
2100 for(RooAbsArg * arg2 : *projDataSel->get()) {
2101 if (!dynamic_cast<RooCategory*>(arg2)) allCat = false ;
2102 }
2103 if (allCat) {
2104 o.binProjData = true ;
2105 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") unbinned projection dataset consist only of discrete variables,"
2106 << " performing projection with binned copy for optimization." << std::endl ;
2107
2108 }
2109 }
2110
2111 // Bin projection dataset if requested
2112 if (o.binProjData) {
2113 projDataSelOwned = std::make_unique<RooDataHist>(std::string(projDataSel->GetName()) + "_binned","Binned projection data",*projDataSel->get(),*projDataSel);
2114 projDataSel = projDataSelOwned.get();
2115 }
2116
2117 // Construct scaled data weighted average
2118 ScaledDataWeightedAverage scaleBind{*projection, *projDataSel, o.scaleFactor, *plotVar};
2119
2120 // Set default range, if not specified
2121 if (o.rangeLo==0 && o.rangeHi==0) {
2122 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2123 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2124 }
2125
2126 // Construct name of curve for data weighed average
2127 std::string curveName(projection->GetName()) ;
2128 curveName.append("_DataAvg[" + projDataSel->get()->contentsString() + "]");
2129 // Append slice set specification if any
2130 if (!sliceSet.empty()) {
2131 curveName.append("_Slice[" + sliceSet.contentsString() + "]");
2132 }
2133 // Append any suffixes imported from RooAbsPdf::plotOn
2134 if (o.curveNameSuffix) {
2135 curveName.append(o.curveNameSuffix) ;
2136 }
2137
2138 // Curve constructor for data weighted average
2140 RooCurve *curve = new RooCurve(projection->GetName(),projection->GetTitle(),scaleBind,
2143
2144 curve->SetName(curveName.c_str()) ;
2145
2146 // Add self to other curve if requested
2147 if (o.addToCurveName) {
2148 RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;
2149
2150 // Curve constructor for sum of curves
2151 RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
2152 sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
2153 delete curve ;
2154 curve = sumCurve ;
2155
2156 }
2157
2158 if (o.curveName) {
2159 curve->SetName(o.curveName) ;
2160 }
2161
2162 // add this new curve to the specified plot frame
2163 frame->addPlotable(curve, o.drawOptions, o.curveInvisible);
2164
2165 } else {
2166
2167 // Set default range, if not specified
2168 if (o.rangeLo==0 && o.rangeHi==0) {
2169 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2170 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2171 }
2172
2173 // Calculate a posteriori range fraction scaling if requested (2nd part of normalization correction for
2174 // result fit on subrange of data)
2175 if (o.postRangeFracScale) {
2176 if (!o.normRangeName) {
2177 o.normRangeName = "plotRange" ;
2178 plotVar->setRange("plotRange",o.rangeLo,o.rangeHi) ;
2179 }
2180
2181 // Evaluate fractional correction integral always on full p.d.f, not component.
2182 GlobalSelectComponentRAII selectCompRAII(true);
2183 std::unique_ptr<RooAbsReal> intFrac{projection->createIntegral(*plotVar,*plotVar,o.normRangeName)};
2184 if(o.stype != RooAbsReal::Raw || this->InheritsFrom(RooAbsPdf::Class())){
2185 // this scaling should only be !=1 when plotting partial ranges
2186 // still, raw means raw
2187 o.scaleFactor /= intFrac->getVal() ;
2188 }
2189 }
2190
2191 // create a new curve of our function using the clone to do the evaluations
2192 // Curve constructor for regular projections
2193
2194 // Set default name of curve
2195 std::string curveName(projection->GetName()) ;
2196 if (!sliceSet.empty()) {
2197 curveName.append("_Slice[" + sliceSet.contentsString() + "]");
2198 }
2199 if (o.curveNameSuffix) {
2200 // Append any suffixes imported from RooAbsPdf::plotOn
2201 curveName.append(o.curveNameSuffix) ;
2202 }
2203
2204 TString opt(o.drawOptions);
2205 if(opt.Contains("P")){
2207 RooHist *graph= new RooHist(*projection,*plotVar,1.,o.scaleFactor,frame->getNormVars(),o.errorFR);
2209
2210 // Override name of curve by user name, if specified
2211 if (o.curveName) {
2212 graph->SetName(o.curveName) ;
2213 }
2214
2215 // add this new curve to the specified plot frame
2217 } else {
2219 RooCurve *curve = new RooCurve(*projection,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
2222 curve->SetName(curveName.c_str()) ;
2223
2224 // Add self to other curve if requested
2225 if (o.addToCurveName) {
2226 RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;
2227 RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
2228 sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
2229 delete curve ;
2230 curve = sumCurve ;
2231 }
2232
2233 // Override name of curve by user name, if specified
2234 if (o.curveName) {
2235 curve->SetName(o.curveName) ;
2236 }
2237
2238 // add this new curve to the specified plot frame
2239 frame->addPlotable(curve, o.drawOptions, o.curveInvisible);
2240 }
2241 }
2242
2243 if (projDataNeededVars) delete projDataNeededVars ;
2244 delete projectionCompList ;
2245 plotVar->setVal(oldPlotVarVal); // reset the plot variable value to not disturb the original state
2246 return frame;
2247}
2248
2249
2250
2251
2252////////////////////////////////////////////////////////////////////////////////
2253/// \deprecated OBSOLETE -- RETAINED FOR BACKWARD COMPATIBILITY. Use plotOn() with Slice() instead
2254
2255RooPlot* RooAbsReal::plotSliceOn(RooPlot *frame, const RooArgSet& sliceSet, Option_t* drawOptions,
2256 double scaleFactor, ScaleType stype, const RooAbsData* projData) const
2257{
2258 RooArgSet projectedVars ;
2259 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
2260
2261 // Take out the sliced variables
2262 for(RooAbsArg * sliceArg : sliceSet) {
2263 RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
2264 if (arg) {
2265 projectedVars.remove(*arg) ;
2266 } else {
2267 coutI(Plotting) << "RooAbsReal::plotSliceOn(" << GetName() << ") slice variable "
2268 << sliceArg->GetName() << " was not projected anyway" << std::endl ;
2269 }
2270 }
2271
2272 PlotOpt o ;
2273 o.drawOptions = drawOptions ;
2274 o.scaleFactor = scaleFactor ;
2275 o.stype = stype ;
2276 o.projData = projData ;
2277 o.projSet = &projectedVars ;
2278 return plotOn(frame,o) ;
2279}
2280
2281
2282
2283
2284//_____________________________________________________________________________
2285// coverity[PASS_BY_VALUE]
2287
2288{
2289 // Plotting engine for asymmetries. Implements the functionality if plotOn(frame,Asymmetry(...)))
2290 //
2291 // Plot asymmetry of ourselves, defined as
2292 //
2293 // asym = f(asymCat=-1) - f(asymCat=+1) / ( f(asymCat=-1) + f(asymCat=+1) )
2294 //
2295 // on frame. If frame contains a histogram, all dimensions of the plotted
2296 // asymmetry function that occur in the previously plotted dataset are projected via partial integration.
2297 // Otherwise no projections are performed,
2298 //
2299 // The asymmetry function can be multiplied with an optional scale factor. The default projection
2300 // behaviour can be overridden by supplying an optional set of dependents to project.
2301
2302 // Sanity checks
2303 if (plotSanityChecks(frame)) return frame ;
2304
2305 // ProjDataVars is either all projData observables, or the user indicated subset of it
2306 RooArgSet projDataVars ;
2307 if (o.projData) {
2308 if (o.projDataSet) {
2309 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(o.projData->get()->selectCommon(*o.projDataSet))};
2310 projDataVars.add(*tmp) ;
2311 } else {
2312 projDataVars.add(*o.projData->get()) ;
2313 }
2314 }
2315
2316 // Must depend on asymCat
2317 if (!dependsOn(asymCat)) {
2318 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2319 << ") function doesn't depend on asymmetry category " << asymCat.GetName() << std::endl ;
2320 return frame ;
2321 }
2322
2323 // asymCat must be a signCat
2324 if (!asymCat.isSignType()) {
2325 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2326 << ") asymmetry category must have 2 or 3 states with index values -1,0,1" << std::endl ;
2327 return frame ;
2328 }
2329
2330 // Make list of variables to be projected
2331 RooArgSet projectedVars ;
2332 RooArgSet sliceSet ;
2333 if (o.projSet) {
2334 makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,false) ;
2335
2336 // Print list of non-projected variables
2337 if (frame->getNormVars()) {
2338 RooArgSet sliceSetTmp;
2339 getObservables(frame->getNormVars(), sliceSetTmp) ;
2340 sliceSetTmp.remove(projectedVars,true,true) ;
2341 sliceSetTmp.remove(*frame->getPlotVar(),true,true) ;
2342
2343 if (o.projData) {
2344 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(projDataVars.selectCommon(*o.projSet))};
2345 sliceSetTmp.remove(*tmp,true,true) ;
2346 }
2347
2348 if (!sliceSetTmp.empty()) {
2349 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on "
2350 << frame->getPlotVar()->GetName() << " represents a slice in " << sliceSetTmp << std::endl ;
2351 }
2352 sliceSet.add(sliceSetTmp) ;
2353 }
2354 } else {
2355 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
2356 }
2357
2358
2359 // Take out data-projected dependens from projectedVars
2360 RooArgSet* projDataNeededVars = nullptr ;
2361 if (o.projData) {
2362 projDataNeededVars = static_cast<RooArgSet*>(projectedVars.selectCommon(projDataVars)) ;
2363 projectedVars.remove(projDataVars,true,true) ;
2364 }
2365
2366 // Take out plotted asymmetry from projection
2367 if (projectedVars.find(asymCat.GetName())) {
2368 projectedVars.remove(*projectedVars.find(asymCat.GetName())) ;
2369 }
2370
2371 // Clone the plot variable
2372 RooAbsReal* realVar = static_cast<RooRealVar*>(frame->getPlotVar()) ;
2373 RooRealVar* plotVar = static_cast<RooRealVar*>(realVar->Clone()) ;
2374
2375 // Inform user about projections
2376 if (!projectedVars.empty()) {
2377 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on " << plotVar->GetName()
2378 << " projects variables " << projectedVars << std::endl ;
2379 }
2380 if (projDataNeededVars && !projDataNeededVars->empty()) {
2381 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2382 << " averages using data variables "<< *projDataNeededVars << std::endl ;
2383 }
2384
2385
2386 // Customize two copies of projection with fixed negative and positive asymmetry
2387 std::unique_ptr<RooAbsCategoryLValue> asymPos{static_cast<RooAbsCategoryLValue*>(asymCat.Clone("asym_pos"))};
2388 std::unique_ptr<RooAbsCategoryLValue> asymNeg{static_cast<RooAbsCategoryLValue*>(asymCat.Clone("asym_neg"))};
2389 asymPos->setIndex(1) ;
2390 asymNeg->setIndex(-1) ;
2391 RooCustomizer custPos{*this,"pos"};
2392 RooCustomizer custNeg{*this,"neg"};
2393 //custPos->setOwning(true) ;
2394 //custNeg->setOwning(true) ;
2395 custPos.replaceArg(asymCat,*asymPos) ;
2396 custNeg.replaceArg(asymCat,*asymNeg) ;
2397 std::unique_ptr<RooAbsReal> funcPos{static_cast<RooAbsReal*>(custPos.build())};
2398 std::unique_ptr<RooAbsReal> funcNeg{static_cast<RooAbsReal*>(custNeg.build())};
2399
2400 // Create projection integral
2401 RooArgSet *posProjCompList;
2402 RooArgSet *negProjCompList;
2403
2404 // Add projDataVars to normalized dependents of projection
2405 // This is needed only for asymmetries (why?)
2406 RooArgSet depPos(*plotVar,*asymPos) ;
2407 RooArgSet depNeg(*plotVar,*asymNeg) ;
2408 depPos.add(projDataVars) ;
2409 depNeg.add(projDataVars) ;
2410
2411 const RooAbsReal *posProj = funcPos->createPlotProjection(depPos, &projectedVars, posProjCompList, o.projectionRangeName) ;
2412 const RooAbsReal *negProj = funcNeg->createPlotProjection(depNeg, &projectedVars, negProjCompList, o.projectionRangeName) ;
2413 if (!posProj || !negProj) {
2414 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") Unable to create projections, abort" << std::endl ;
2415 return frame ;
2416 }
2417
2418 // Create a RooFormulaVar representing the asymmetry
2419 TString asymName(GetName()) ;
2420 asymName.Append("_Asym[") ;
2421 asymName.Append(asymCat.GetName()) ;
2422 asymName.Append("]") ;
2423 TString asymTitle(asymCat.GetName()) ;
2424 asymTitle.Append(" Asymmetry of ") ;
2425 asymTitle.Append(GetTitle()) ;
2426 RooFormulaVar funcAsym{asymName,asymTitle,"(@0-@1)/(@0+@1)",RooArgSet(*posProj,*negProj)};
2427
2428 if (o.projData) {
2429
2430 // If data set contains more rows than needed, make reduced copy first
2431 RooAbsData* projDataSel = const_cast<RooAbsData*>(o.projData);
2432 std::unique_ptr<RooAbsData> projDataSelOwned;
2433 if (projDataNeededVars && projDataNeededVars->size() < o.projData->get()->size()) {
2434
2435 // Determine if there are any slice variables in the projection set
2436 RooArgSet sliceDataSet;
2437 sliceSet.selectCommon(*o.projData->get(), sliceDataSet);
2438 TString cutString ;
2439 if (!sliceDataSet.empty()) {
2440 bool first(true) ;
2441 for(RooAbsArg * sliceVar : sliceDataSet) {
2442 if (!first) {
2443 cutString.Append("&&") ;
2444 } else {
2445 first=false ;
2446 }
2447
2448 RooAbsRealLValue* real ;
2450 if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
2451 cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
2452 } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
2453 cutString.Append(Form("%s==%d",cat->GetName(),cat->getCurrentIndex())) ;
2454 }
2455 }
2456 }
2457
2458 if (!cutString.IsNull()) {
2459 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2460 << ") reducing given projection dataset to entries with " << cutString << std::endl ;
2461 }
2462 projDataSelOwned = std::unique_ptr<RooAbsData>{const_cast<RooAbsData*>(o.projData)->reduce(*projDataNeededVars,cutString.IsNull() ? nullptr : cutString)};
2463 projDataSel = projDataSelOwned.get();
2464 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2465 << ") only the following components of the projection data will be used: " << *projDataNeededVars << std::endl ;
2466 }
2467
2468
2469 // Construct scaled data weighted average
2470 ScaledDataWeightedAverage scaleBind{funcAsym, *projDataSel, o.scaleFactor, *plotVar};
2471
2472 // Set default range, if not specified
2473 if (o.rangeLo==0 && o.rangeHi==0) {
2474 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2475 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2476 }
2477
2478 // Construct name of curve for data weighed average
2479 TString curveName(funcAsym.GetName()) ;
2480 curveName.Append(Form("_DataAvg[%s]",projDataSel->get()->contentsString().c_str())) ;
2481 // Append slice set specification if any
2482 if (!sliceSet.empty()) {
2483 curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
2484 }
2485 // Append any suffixes imported from RooAbsPdf::plotOn
2486 if (o.curveNameSuffix) {
2487 curveName.Append(o.curveNameSuffix) ;
2488 }
2489
2490
2492 RooCurve *curve = new RooCurve(funcAsym.GetName(),funcAsym.GetTitle(),scaleBind,
2493 o.rangeLo,o.rangeHi,frame->GetNbinsX(),o.precision,o.precision,false,o.wmode,o.numee,o.doeeval,o.eeval) ;
2495
2496 dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;
2497 // add this new curve to the specified plot frame
2498 frame->addPlotable(curve, o.drawOptions);
2499
2500 ccoutW(Eval) << std::endl ;
2501 } else {
2502
2503 // Set default range, if not specified
2504 if (o.rangeLo==0 && o.rangeHi==0) {
2505 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2506 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2507 }
2508
2510 RooCurve* curve= new RooCurve(funcAsym,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
2511 o.scaleFactor,nullptr,o.precision,o.precision,false,o.wmode,o.numee,o.doeeval,o.eeval);
2513
2514 dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;
2515
2516
2517 // Set default name of curve
2518 TString curveName(funcAsym.GetName()) ;
2519 if (!sliceSet.empty()) {
2520 curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
2521 }
2522 if (o.curveNameSuffix) {
2523 // Append any suffixes imported from RooAbsPdf::plotOn
2524 curveName.Append(o.curveNameSuffix) ;
2525 }
2526 curve->SetName(curveName.Data()) ;
2527
2528 // add this new curve to the specified plot frame
2529 frame->addPlotable(curve, o.drawOptions);
2530
2531 }
2532
2533 // Cleanup
2534 delete posProjCompList ;
2535 delete negProjCompList ;
2536
2537 delete plotVar ;
2538
2539 return frame;
2540}
2541
2542
2543
2544////////////////////////////////////////////////////////////////////////////////
2545/// \brief Propagates parameter uncertainties to an uncertainty estimate for this RooAbsReal.
2546///
2547/// Estimates the uncertainty \f$\sigma_f(x;\theta)\f$ on a function \f$f(x;\theta)\f$ represented by this RooAbsReal.
2548/// Here, \f$\theta\f$ is a vector of parameters with uncertainties \f$\sigma_\theta\f$, and \f$x\f$ are usually observables.
2549/// The uncertainty is estimated by *linearly* propagating the parameter uncertainties using the correlation matrix from a fit result.
2550///
2551/// The square of the uncertainty on \f$f(x;\theta)\f$ is calculated as follows:
2552/// \f[
2553/// \sigma_f(x)^2 = \Delta f_i(x) \cdot \mathrm{Corr}_{i, j} \cdot \Delta f_j(x),
2554/// \f]
2555/// where \f$ \Delta f_i(x) = \frac{1}{2} \left(f(x;\theta_i + \sigma_{\theta_i}) - f(x; \theta_i - \sigma_{\theta_i}) \right) \f$
2556/// is the vector of function variations when changing the parameters one at a time, and
2557/// \f$ \mathrm{Corr}_{i,j} = \left(\sigma_{\theta_i} \sigma_{\theta_j}\right)^{-1} \cdot \mathrm{Cov}_{i,j} \f$ is the correlation matrix from the fit result.
2558
2559double RooAbsReal::getPropagatedError(const RooFitResult &fr, const RooArgSet &nset) const
2560{
2561 // Calling getParameters() might be costly, but necessary to get the right
2562 // parameters in the RooAbsReal. The RooFitResult only stores snapshots.
2563 RooArgSet allParamsInAbsReal;
2564 getParameters(&nset, allParamsInAbsReal);
2565
2566 RooArgList paramList;
2567 for(auto * rrvFitRes : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {
2568
2569 auto rrvInAbsReal = static_cast<RooRealVar const*>(allParamsInAbsReal.find(*rrvFitRes));
2570
2571 // If this RooAbsReal is a RooRealVar in the fit result, we don't need to
2572 // propagate anything and can just return the error in the fit result
2573 if(rrvFitRes->namePtr() == namePtr()) return rrvFitRes->getError();
2574
2575 // Strip out parameters with zero error
2576 if (rrvFitRes->getError() <= std::abs(rrvFitRes->getVal()) * std::numeric_limits<double>::epsilon()) continue;
2577
2578 // Ignore parameters in the fit result that this RooAbsReal doesn't depend on
2579 if(!rrvInAbsReal) continue;
2580
2581 // Checking for float equality is a bad. We check if the values are
2582 // negligibly far away from each other, relative to the uncertainty.
2583 if(std::abs(rrvInAbsReal->getVal() - rrvFitRes->getVal()) > 0.01 * rrvFitRes->getError()) {
2584 std::stringstream errMsg;
2585 errMsg << "RooAbsReal::getPropagatedError(): the parameters of the RooAbsReal don't have"
2586 << " the same values as in the fit result! The logic of getPropagatedError is broken in this case.";
2587
2588 throw std::runtime_error(errMsg.str());
2589 }
2590
2591 paramList.add(*rrvInAbsReal);
2592 }
2593
2594 std::vector<double> plusVar;
2595 std::vector<double> minusVar;
2596 plusVar.reserve(paramList.size());
2597 minusVar.reserve(paramList.size());
2598
2599 // Create std::vector of plus,minus variations for each parameter
2600 TMatrixDSym V(paramList.size() == fr.floatParsFinal().size() ?
2601 fr.covarianceMatrix() :
2602 fr.reducedCovarianceMatrix(paramList)) ;
2603
2604 for (std::size_t ivar=0 ; ivar<paramList.size() ; ivar++) {
2605
2606 auto& rrv = static_cast<RooRealVar&>(paramList[ivar]);
2607
2608 double cenVal = rrv.getVal() ;
2609 double errVal = sqrt(V(ivar,ivar)) ;
2610
2611 // Make Plus variation
2612 rrv.setVal(cenVal+errVal) ;
2613 plusVar.push_back(getVal(nset)) ;
2614
2615 // Make Minus variation
2616 rrv.setVal(cenVal-errVal) ;
2617 minusVar.push_back(getVal(nset)) ;
2618
2619 rrv.setVal(cenVal) ;
2620 }
2621
2622 // Re-evaluate this RooAbsReal with the central parameters just to be
2623 // extra-safe that a call to `getPropagatedError()` doesn't change any state.
2624 // It should not be necessary because thanks to the dirty flag propagation
2625 // the RooAbsReal is re-evaluated anyway the next time getVal() is called.
2626 // Still there are imaginable corner cases where it would not be triggered,
2627 // for example if the user changes the RooFit operation more after the error
2628 // propagation.
2629 getVal(nset);
2630
2631 TMatrixDSym C(paramList.size()) ;
2632 std::vector<double> errVec(paramList.size()) ;
2633 for (std::size_t i=0 ; i<paramList.size() ; i++) {
2634 errVec[i] = std::sqrt(V(i,i)) ;
2635 for (std::size_t j=i ; j<paramList.size() ; j++) {
2636 C(i,j) = V(i,j) / std::sqrt(V(i,i)*V(j,j));
2637 C(j,i) = C(i,j) ;
2638 }
2639 }
2640
2641 // Make std::vector of variations
2642 TVectorD F(plusVar.size()) ;
2643 for (std::size_t j=0 ; j<plusVar.size() ; j++) {
2644 F[j] = (plusVar[j]-minusVar[j]) * 0.5;
2645 }
2646
2647 // Calculate error in linear approximation from variations and correlation coefficient
2648 double sum = F*(C*F) ;
2649
2650 return sqrt(sum) ;
2651}
2652
2653
2654
2655////////////////////////////////////////////////////////////////////////////////
2656/// Plot function or PDF on frame with support for visualization of the uncertainty encoded in the given fit result fr.
2657/// \param[in] frame RooPlot to plot on
2658/// \param[in] fr The RooFitResult, where errors can be extracted
2659/// \param[in] Z The desired significance (width) of the error band
2660/// \param[in] params If non-zero, consider only the subset of the parameters in fr for the error evaluation
2661/// \param[in] argList Optional `RooCmdArg` that can be applied to a regular plotOn() operation
2662/// \param[in] linMethod By default (linMethod=true), a linearized error is shown.
2663/// \return The RooPlot the band was plotted on (for chaining of plotting commands).
2664///
2665/// The linearized error is calculated as follows:
2666/// \f[
2667/// \mathrm{error}(x) = Z * F_a(x) * \mathrm{Corr}(a,a') * F_{a'}^\mathrm{T}(x),
2668/// \f]
2669///
2670/// where
2671/// \f[
2672/// F_a(x) = \frac{ f(x,a+\mathrm{d}a) - f(x,a-\mathrm{d}a) }{2},
2673/// \f]
2674/// with \f$ f(x) \f$ the plotted curve and \f$ \mathrm{d}a \f$ taken from the fit result, and
2675/// \f$ \mathrm{Corr}(a,a') \f$ = the correlation matrix from the fit result, and \f$ Z \f$ = requested signifance (\f$ Z \sigma \f$ band)
2676///
2677/// The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters), but may
2678/// not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made
2679///
2680/// Alternatively, a more robust error is calculated using a sampling method. In this method a number of curves
2681/// is calculated with variations of the parameter values, as drawn from a multi-variate Gaussian p.d.f. that is constructed
2682/// from the fit results covariance matrix. The error(x) is determined by calculating a central interval that capture N% of the variations
2683/// for each value 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
2684/// 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)
2685/// Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations, but may take (much)
2686/// longer to calculate.
2687
2688RooPlot* RooAbsReal::plotOnWithErrorBand(RooPlot* frame,const RooFitResult& fr, double Z,const RooArgSet* params, const RooLinkedList& argList, bool linMethod) const
2689{
2690 RooLinkedList plotArgListTmp(argList) ;
2691 RooCmdConfig::stripCmdList(plotArgListTmp,"VisualizeError,MoveToBack") ;
2692
2693 // Strip any 'internal normalization' arguments from list
2694 RooLinkedList plotArgList ;
2695 for (auto * cmd : static_range_cast<RooCmdArg*>(plotArgListTmp)) {
2696 if (std::string("Normalization")==cmd->GetName()) {
2697 if (((RooCmdArg*)cmd)->getInt(1)!=0) {
2698 } else {
2699 plotArgList.Add(cmd) ;
2700 }
2701 } else {
2702 plotArgList.Add(cmd) ;
2703 }
2704 }
2705
2706 // Function to plot a single curve, creating a copy of the plotArgList to
2707 // pass as plot command arguments. The "FillColor" command is removed because
2708 // it has no effect on plotting single curves and would cause a warning.
2709 auto plotFunc = [&](RooAbsReal const& absReal) {
2710 RooLinkedList tmp(plotArgList) ;
2711 RooCmdConfig::stripCmdList(tmp, "FillColor");
2712 absReal.plotOn(frame, tmp);
2713 };
2714
2715 // Generate central value curve
2716 plotFunc(*this);
2717 RooCurve* cenCurve = frame->getCurve() ;
2718 if(!cenCurve){
2719 coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOnWithErrorBand: no curve for central value available" << std::endl;
2720 return frame;
2721 }
2722 frame->remove(nullptr,false) ;
2723
2724 RooCurve* band(nullptr) ;
2725 if (!linMethod) {
2726
2727 // *** Interval method ***
2728 //
2729 // Make N variations of parameters samples from V and visualize N% central interval where N% is defined from Z
2730
2731 // Clone self for internal use
2732 RooAbsReal* cloneFunc = static_cast<RooAbsReal*>(cloneTree()) ;
2733 RooArgSet cloneParams;
2734 cloneFunc->getObservables(&fr.floatParsFinal(), cloneParams) ;
2735 RooArgSet errorParams{cloneParams};
2736 if(params) {
2737 // clear and fill errorParams only with parameters that both in params and cloneParams
2738 cloneParams.selectCommon(*params, errorParams);
2739 }
2740
2741 // Generate 100 random parameter points distributed according to fit result covariance matrix
2742 RooAbsPdf* paramPdf = fr.createHessePdf(errorParams) ;
2743 Int_t n = Int_t(100./TMath::Erfc(Z/sqrt(2.))) ;
2744 if (n<100) n=100 ;
2745
2746 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") INFO: visualizing " << Z << "-sigma uncertainties in parameters "
2747 << errorParams << " from fit result " << fr.GetName() << " using " << n << " samplings." << std::endl ;
2748
2749 // Generate variation curves with above set of parameter values
2750 double ymin = frame->GetMinimum() ;
2751 double ymax = frame->GetMaximum() ;
2752 std::unique_ptr<RooDataSet> generatedData{paramPdf->generate(errorParams,n)};
2753 std::vector<RooCurve*> cvec ;
2754 for (int i=0 ; i<generatedData->numEntries() ; i++) {
2755 cloneParams.assign(*generatedData->get(i)) ;
2756 plotFunc(*cloneFunc);
2757 cvec.push_back(frame->getCurve()) ;
2758 frame->remove(nullptr,false) ;
2759 }
2760 frame->SetMinimum(ymin) ;
2761 frame->SetMaximum(ymax) ;
2762
2763
2764 // Generate upper and lower curve points from 68% interval around each point of central curve
2765 band = cenCurve->makeErrorBand(cvec,Z) ;
2766
2767 // Cleanup
2768 delete paramPdf ;
2769 delete cloneFunc ;
2770 for (std::vector<RooCurve*>::iterator i=cvec.begin() ; i!=cvec.end() ; ++i) {
2771 delete (*i) ;
2772 }
2773
2774 } else {
2775
2776 // *** Linear Method ***
2777 //
2778 // Make a one-sigma up- and down fluctation for each parameter and visualize
2779 // a from a linearized calculation as follows
2780 //
2781 // error(x) = F(a) C_aa' F(a')
2782 //
2783 // Where F(a) = (f(x,a+da) - f(x,a-da))/2
2784 // and C_aa' is the correlation matrix
2785
2786 // Strip out parameters with zero error
2787 RooArgList fpf_stripped;
2788 for (auto const* frv : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {
2789 if (frv->getError() > frv->getVal() * std::numeric_limits<double>::epsilon()) {
2790 fpf_stripped.add(*frv);
2791 }
2792 }
2793
2794 // Clone self for internal use
2795 RooAbsReal* cloneFunc = static_cast<RooAbsReal*>(cloneTree()) ;
2796 RooArgSet cloneParams;
2797 cloneFunc->getObservables(&fpf_stripped, cloneParams) ;
2798 RooArgSet errorParams{cloneParams};
2799 if(params) {
2800 // clear and fill errorParams only with parameters that both in params and cloneParams
2801 cloneParams.selectCommon(*params, errorParams);
2802 }
2803
2804
2805 // Make list of parameter instances of cloneFunc in order of error matrix
2806 RooArgList paramList ;
2807 const RooArgList& fpf = fr.floatParsFinal() ;
2808 std::vector<int> fpf_idx ;
2809 for (std::size_t i=0 ; i<fpf.size() ; i++) {
2810 RooAbsArg* par = errorParams.find(fpf[i].GetName()) ;
2811 if (par) {
2812 paramList.add(*par) ;
2813 fpf_idx.push_back(i) ;
2814 }
2815 }
2816
2817 std::vector<RooCurve *> plusVar;
2818 std::vector<RooCurve *> minusVar;
2819
2820 // Create std::vector of plus,minus variations for each parameter
2821
2822 TMatrixDSym V(paramList.size() == fr.floatParsFinal().size() ?
2823 fr.covarianceMatrix():
2824 fr.reducedCovarianceMatrix(paramList)) ;
2825
2826
2827 for (std::size_t ivar=0 ; ivar<paramList.size() ; ivar++) {
2828
2829 RooRealVar& rrv = static_cast<RooRealVar&>(fpf[fpf_idx[ivar]]) ;
2830
2831 double cenVal = rrv.getVal() ;
2832 double errVal = sqrt(V(ivar,ivar)) ;
2833
2834 // Make Plus variation
2835 (static_cast<RooRealVar*>(paramList.at(ivar)))->setVal(cenVal+Z*errVal) ;
2836
2837
2838 plotFunc(*cloneFunc);
2839 plusVar.push_back(frame->getCurve()) ;
2840 frame->remove(nullptr,false) ;
2841
2842
2843 // Make Minus variation
2844 (static_cast<RooRealVar*>(paramList.at(ivar)))->setVal(cenVal-Z*errVal) ;
2845 plotFunc(*cloneFunc);
2846 minusVar.push_back(frame->getCurve()) ;
2847 frame->remove(nullptr,false) ;
2848
2849 (static_cast<RooRealVar*>(paramList.at(ivar)))->setVal(cenVal) ;
2850 }
2851
2852 TMatrixDSym C(paramList.size()) ;
2853 std::vector<double> errVec(paramList.size()) ;
2854 for (std::size_t i=0 ; i<paramList.size() ; i++) {
2855 errVec[i] = sqrt(V(i,i)) ;
2856 for (std::size_t j=i ; j<paramList.size() ; j++) {
2857 C(i,j) = V(i,j)/sqrt(V(i,i)*V(j,j)) ;
2858 C(j,i) = C(i,j) ;
2859 }
2860 }
2861
2862 band = cenCurve->makeErrorBand(plusVar,minusVar,C,Z) ;
2863
2864
2865 // Cleanup
2866 delete cloneFunc ;
2867 for (std::vector<RooCurve*>::iterator i=plusVar.begin() ; i!=plusVar.end() ; ++i) {
2868 delete (*i) ;
2869 }
2870 for (std::vector<RooCurve*>::iterator i=minusVar.begin() ; i!=minusVar.end() ; ++i) {
2871 delete (*i) ;
2872 }
2873
2874 }
2875
2876 delete cenCurve ;
2877 if (!band) return frame ;
2878
2879 // Define configuration for this method
2880 RooCmdConfig pc("RooAbsPdf::plotOn(" + std::string(GetName()) + ")");
2881 pc.defineString("drawOption","DrawOption",0,"F") ;
2882 pc.defineString("curveNameSuffix","CurveNameSuffix",0,"") ;
2883 pc.defineInt("lineColor","LineColor",0,-999) ;
2884 pc.defineInt("lineStyle","LineStyle",0,-999) ;
2885 pc.defineInt("lineWidth","LineWidth",0,-999) ;
2886 pc.defineInt("markerColor","MarkerColor",0,-999) ;
2887 pc.defineInt("markerStyle","MarkerStyle",0,-999) ;
2888 pc.defineDouble("markerSize","MarkerSize",0,-999) ;
2889 pc.defineInt("fillColor","FillColor",0,-999) ;
2890 pc.defineInt("fillStyle","FillStyle",0,-999) ;
2891 pc.defineString("curveName","Name",0,"") ;
2892 pc.defineInt("curveInvisible","Invisible",0,0) ;
2893 pc.defineInt("moveToBack","MoveToBack",0,0) ;
2894 pc.allowUndefined() ;
2895
2896 // Process & check varargs
2897 pc.process(argList) ;
2898 if (!pc.ok(true)) {
2899 return frame ;
2900 }
2901
2902 // Insert error band in plot frame
2903 frame->addPlotable(band,pc.getString("drawOption"),pc.getInt("curveInvisible")) ;
2904
2905 // Optionally adjust line/fill attributes
2906 Int_t lineColor = pc.getInt("lineColor") ;
2907 Int_t lineStyle = pc.getInt("lineStyle") ;
2908 Int_t lineWidth = pc.getInt("lineWidth") ;
2909 Int_t markerColor = pc.getInt("markerColor") ;
2910 Int_t markerStyle = pc.getInt("markerStyle") ;
2911 Size_t markerSize = pc.getDouble("markerSize") ;
2912 Int_t fillColor = pc.getInt("fillColor") ;
2913 Int_t fillStyle = pc.getInt("fillStyle") ;
2914 if (lineColor!=-999) frame->getAttLine()->SetLineColor(lineColor) ;
2915 if (lineStyle!=-999) frame->getAttLine()->SetLineStyle(lineStyle) ;
2916 if (lineWidth!=-999) frame->getAttLine()->SetLineWidth(lineWidth) ;
2917 if (fillColor!=-999) frame->getAttFill()->SetFillColor(fillColor) ;
2918 if (fillStyle!=-999) frame->getAttFill()->SetFillStyle(fillStyle) ;
2919 if (markerColor!=-999) frame->getAttMarker()->SetMarkerColor(markerColor) ;
2920 if (markerStyle!=-999) frame->getAttMarker()->SetMarkerStyle(markerStyle) ;
2921 if (markerSize!=-999) frame->getAttMarker()->SetMarkerSize(markerSize) ;
2922
2923 // Adjust name if requested
2924 if (pc.getString("curveName",nullptr,true)) {
2925 band->SetName(pc.getString("curveName",nullptr,true)) ;
2926 } else if (pc.getString("curveNameSuffix",nullptr,true)) {
2927 TString name(band->GetName()) ;
2928 name.Append(pc.getString("curveNameSuffix",nullptr,true)) ;
2929 band->SetName(name.Data()) ;
2930 }
2931
2932 // Move last inserted object to back to drawing stack if requested
2933 if (pc.getInt("moveToBack") && frame->numItems()>1) {
2934 frame->drawBefore(frame->getObject(0)->GetName(), frame->getCurve()->GetName());
2935 }
2936
2937
2938 return frame ;
2939}
2940
2941
2942
2943
2944////////////////////////////////////////////////////////////////////////////////
2945/// Utility function for plotOn(), perform general sanity check on frame to ensure safe plotting operations
2946
2948{
2949 // check that we are passed a valid plot frame to use
2950 if(nullptr == frame) {
2951 coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << std::endl;
2952 return true;
2953 }
2954
2955 // check that this frame knows what variable to plot
2956 RooAbsReal* var = frame->getPlotVar() ;
2957 if(!var) {
2958 coutE(Plotting) << ClassName() << "::" << GetName()
2959 << ":plotOn: frame does not specify a plot variable" << std::endl;
2960 return true;
2961 }
2962
2963 // check that the plot variable is not derived
2964 if(!dynamic_cast<RooAbsRealLValue*>(var)) {
2965 coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: cannot plot variable \""
2966 << var->GetName() << "\" of type " << var->ClassName() << std::endl;
2967 return true;
2968 }
2969
2970 // check if we actually depend on the plot variable
2971 if(!this->dependsOn(*var)) {
2972 coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOn: WARNING: variable is not an explicit dependent: "
2973 << var->GetName() << std::endl;
2974 }
2975
2976 return false ;
2977}
2978
2979
2980
2981
2982////////////////////////////////////////////////////////////////////////////////
2983/// Utility function for plotOn() that constructs the set of
2984/// observables to project when plotting ourselves as function of
2985/// 'plotVar'. 'allVars' is the list of variables that must be
2986/// projected, but may contain variables that we do not depend on. If
2987/// 'silent' is cleared, warnings about inconsistent input parameters
2988/// will be printed.
2989
2990void RooAbsReal::makeProjectionSet(const RooAbsArg* plotVar, const RooArgSet* allVars,
2991 RooArgSet& projectedVars, bool silent) const
2992{
2993 cxcoutD(Plotting) << "RooAbsReal::makeProjectionSet(" << GetName() << ") plotVar = " << plotVar->GetName()
2994 << " allVars = " << (allVars?(*allVars):RooArgSet()) << std::endl ;
2995
2996 projectedVars.removeAll() ;
2997 if (!allVars) return ;
2998
2999 // Start out with suggested list of variables
3000 projectedVars.add(*allVars) ;
3001
3002 // Take out plot variable
3003 RooAbsArg *found= projectedVars.find(plotVar->GetName());
3004 if(found) {
3005 projectedVars.remove(*found);
3006
3007 // Take out eventual servers of plotVar
3008 std::unique_ptr<RooArgSet> plotServers{plotVar->getObservables(&projectedVars)};
3009 for(RooAbsArg * ps : *plotServers) {
3010 RooAbsArg* tmp = projectedVars.find(ps->GetName()) ;
3011 if (tmp) {
3012 cxcoutD(Plotting) << "RooAbsReal::makeProjectionSet(" << GetName() << ") removing " << tmp->GetName()
3013 << " from projection set because it a server of " << plotVar->GetName() << std::endl ;
3014 projectedVars.remove(*tmp) ;
3015 }
3016 }
3017
3018 if (!silent) {
3019 coutW(Plotting) << "RooAbsReal::plotOn(" << GetName()
3020 << ") WARNING: cannot project out frame variable ("
3021 << found->GetName() << "), ignoring" << std::endl ;
3022 }
3023 }
3024
3025 // Take out all non-dependents of function
3026 for(RooAbsArg * arg : *allVars) {
3027 if (!dependsOnValue(*arg)) {
3028 projectedVars.remove(*arg,true) ;
3029
3030 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName()
3031 << ") function doesn't depend on projection variable "
3032 << arg->GetName() << ", ignoring" << std::endl ;
3033 }
3034 }
3035}
3036
3037
3038
3039
3040////////////////////////////////////////////////////////////////////////////////
3041/// If true, the current pdf is a selected component (for use in plotting)
3042
3044{
3045 return _selectComp || _globalSelectComp ;
3046}
3047
3048
3049
3050////////////////////////////////////////////////////////////////////////////////
3051/// Global switch controlling the activation of the selectComp() functionality
3052
3054{
3055 _globalSelectComp = flag ;
3056}
3057
3058
3059
3060
3061////////////////////////////////////////////////////////////////////////////////
3062/// Create an interface adaptor f(vars) that binds us to the specified variables
3063/// (in arbitrary order). For example, calling bindVars({x1,x3}) on an object
3064/// F(x1,x2,x3,x4) returns an object f(x1,x3) that is evaluated using the
3065/// current values of x2 and x4. The caller takes ownership of the returned adaptor.
3066
3067RooFit::OwningPtr<RooAbsFunc> RooAbsReal::bindVars(const RooArgSet &vars, const RooArgSet* nset, bool clipInvalid) const
3068{
3069 auto binding = std::make_unique<RooRealBinding>(*this,vars,nset,clipInvalid);
3070 if(!binding->isValid()) {
3071 coutE(InputArguments) << ClassName() << "::" << GetName() << ":bindVars: cannot bind to " << vars << std::endl ;
3072 return nullptr;
3073 }
3074 return RooFit::makeOwningPtr(std::unique_ptr<RooAbsFunc>{std::move(binding)});
3075}
3076
3077
3078
3079////////////////////////////////////////////////////////////////////////////////
3080/// Copy the cached value of another RooAbsArg to our cache.
3081/// Warning: This function just copies the cached values of source,
3082/// it is the callers responsibility to make sure the cache is clean.
3083
3084void RooAbsReal::copyCache(const RooAbsArg* source, bool /*valueOnly*/, bool setValDirty)
3085{
3086 auto other = static_cast<const RooAbsReal*>(source);
3087 assert(dynamic_cast<const RooAbsReal*>(source));
3088
3089 _value = other->_treeReadBuffer ? other->_treeReadBuffer->operator double() : other->_value;
3090
3091 if (setValDirty) {
3092 setValueDirty() ;
3093 }
3094}
3095
3096
3097////////////////////////////////////////////////////////////////////////////////
3098
3100{
3101 vstore.addReal(this)->setBuffer(this,&_value);
3102}
3103
3104
3105////////////////////////////////////////////////////////////////////////////////
3106/// Attach object to a branch of given TTree. By default it will
3107/// register the internal value cache RooAbsReal::_value as branch
3108/// buffer for a double tree branch with the same name as this
3109/// object. If no double branch is found with the name of this
3110/// object, this method looks for a Float_t Int_t, UChar_t and UInt_t, etc
3111/// branch. If any of these are found, a TreeReadBuffer
3112/// that branch is created, and saved in _treeReadBuffer.
3113/// TreeReadBuffer::operator double() can be used to convert the values.
3114/// This is used by copyCache().
3116{
3117 // First determine if branch is taken
3118 TString cleanName(cleanBranchName()) ;
3119 TBranch* branch = t.GetBranch(cleanName) ;
3120 if (branch) {
3121
3122 // Determine if existing branch is Float_t or double
3123 TLeaf* leaf = static_cast<TLeaf*>(branch->GetListOfLeaves()->At(0)) ;
3124
3125 // Check that leaf is _not_ an array
3126 Int_t dummy ;
3127 TLeaf* counterLeaf = leaf->GetLeafCounter(dummy) ;
3128 if (counterLeaf) {
3129 coutE(Eval) << "RooAbsReal::attachToTree(" << GetName() << ") ERROR: TTree branch " << GetName()
3130 << " is an array and cannot be attached to a RooAbsReal" << std::endl ;
3131 return ;
3132 }
3133
3134 TString typeName(leaf->GetTypeName()) ;
3135
3136
3137 // For different type names, store three items:
3138 // first: A tag attached to this instance. Not used inside RooFit, any more, but users might rely on it.
3139 // second: A function to attach
3140 std::map<std::string, std::pair<std::string, std::function<std::unique_ptr<TreeReadBuffer>()>>> typeMap {
3141 {"Float_t", {"FLOAT_TREE_BRANCH", [&](){ return createTreeReadBuffer<Float_t >(cleanName, t); }}},
3142 {"Int_t", {"INTEGER_TREE_BRANCH", [&](){ return createTreeReadBuffer<Int_t >(cleanName, t); }}},
3143 {"UChar_t", {"BYTE_TREE_BRANCH", [&](){ return createTreeReadBuffer<UChar_t >(cleanName, t); }}},
3144 {"Bool_t", {"BOOL_TREE_BRANCH", [&](){ return createTreeReadBuffer<Bool_t >(cleanName, t); }}},
3145 {"Char_t", {"SIGNEDBYTE_TREE_BRANCH", [&](){ return createTreeReadBuffer<Char_t >(cleanName, t); }}},
3146 {"UInt_t", {"UNSIGNED_INTEGER_TREE_BRANCH", [&](){ return createTreeReadBuffer<UInt_t >(cleanName, t); }}},
3147 {"Long64_t", {"LONG_TREE_BRANCH", [&](){ return createTreeReadBuffer<Long64_t >(cleanName, t); }}},
3148 {"ULong64_t", {"UNSIGNED_LONG_TREE_BRANCH", [&](){ return createTreeReadBuffer<ULong64_t>(cleanName, t); }}},
3149 {"Short_t", {"SHORT_TREE_BRANCH", [&](){ return createTreeReadBuffer<Short_t >(cleanName, t); }}},
3150 {"UShort_t", {"UNSIGNED_SHORT_TREE_BRANCH", [&](){ return createTreeReadBuffer<UShort_t >(cleanName, t); }}},
3151 };
3152
3153 auto typeDetails = typeMap.find(typeName.Data());
3154 if (typeDetails != typeMap.end()) {
3155 coutI(DataHandling) << "RooAbsReal::attachToTree(" << GetName() << ") TTree " << typeDetails->first << " branch " << GetName()
3156 << " will be converted to double precision." << std::endl ;
3157 setAttribute(typeDetails->second.first.c_str(), true);
3158 _treeReadBuffer = typeDetails->second.second();
3159 } else {
3160 _treeReadBuffer = nullptr;
3161
3162 if (!typeName.CompareTo("Double_t")) {
3163 t.SetBranchAddress(cleanName, &_value);
3164 }
3165 else {
3166 coutE(InputArguments) << "RooAbsReal::attachToTree(" << GetName() << ") data type " << typeName << " is not supported." << std::endl ;
3167 }
3168 }
3169 } else {
3170
3171 TString format(cleanName);
3172 format.Append("/D");
3173 branch = t.Branch(cleanName, &_value, (const Text_t*)format, bufSize);
3174 }
3175
3176}
3177
3178
3179
3180////////////////////////////////////////////////////////////////////////////////
3181/// Fill the tree branch that associated with this object with its current value
3182
3184{
3185 // First determine if branch is taken
3186 TBranch* branch = t.GetBranch(cleanBranchName()) ;
3187 if (!branch) {
3188 coutE(Eval) << "RooAbsReal::fillTreeBranch(" << GetName() << ") ERROR: not attached to tree: " << cleanBranchName() << std::endl ;
3189 assert(0) ;
3190 }
3191 branch->Fill() ;
3192
3193}
3194
3195
3196
3197////////////////////////////////////////////////////////////////////////////////
3198/// (De)Activate associated tree branch
3199
3201{
3202 TBranch* branch = t.GetBranch(cleanBranchName()) ;
3203 if (branch) {
3204 t.SetBranchStatus(cleanBranchName(),active?true:false) ;
3205 }
3206}
3207
3208
3209
3210////////////////////////////////////////////////////////////////////////////////
3211/// Create a RooRealVar fundamental object with our properties. The new
3212/// object will be created without any fit limits.
3213
3215{
3216 auto fund = std::make_unique<RooRealVar>(newname?newname:GetName(),GetTitle(),_value,getUnit());
3217 fund->removeRange();
3218 fund->setPlotLabel(getPlotLabel());
3219 fund->setAttribute("fundamentalCopy");
3220 return RooFit::makeOwningPtr<RooAbsArg>(std::move(fund));
3221}
3222
3223
3224
3225////////////////////////////////////////////////////////////////////////////////
3226/// Utility function for use in getAnalyticalIntegral(). If the
3227/// content of proxy 'a' occurs in set 'allDeps' then the argument
3228/// held in 'a' is copied from allDeps to analDeps
3229
3230bool RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
3231 const RooArgProxy& a) const
3232{
3233 TList nameList ;
3234 nameList.Add(new TObjString(a.absArg()->GetName())) ;
3235 bool result = matchArgsByName(allDeps,analDeps,nameList) ;
3236 nameList.Delete() ;
3237 return result ;
3238}
3239
3240
3241
3242////////////////////////////////////////////////////////////////////////////////
3243/// Utility function for use in getAnalyticalIntegral(). If the
3244/// contents of proxies a,b occur in set 'allDeps' then the arguments
3245/// held in a,b are copied from allDeps to analDeps
3246
3247bool RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
3248 const RooArgProxy& a, const RooArgProxy& b) const
3249{
3250 TList nameList ;
3251 nameList.Add(new TObjString(a.absArg()->GetName())) ;
3252 nameList.Add(new TObjString(b.absArg()->GetName())) ;
3253 bool result = matchArgsByName(allDeps,analDeps,nameList) ;
3254 nameList.Delete() ;
3255 return result ;
3256}
3257
3258
3259
3260////////////////////////////////////////////////////////////////////////////////
3261/// Utility function for use in getAnalyticalIntegral(). If the
3262/// contents of proxies a,b,c occur in set 'allDeps' then the arguments
3263/// held in a,b,c are copied from allDeps to analDeps
3264
3265bool RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
3266 const RooArgProxy& a, const RooArgProxy& b,
3267 const RooArgProxy& c) const
3268{
3269 TList nameList ;
3270 nameList.Add(new TObjString(a.absArg()->GetName())) ;
3271 nameList.Add(new TObjString(b.absArg()->GetName())) ;
3272 nameList.Add(new TObjString(c.absArg()->GetName())) ;
3273 bool result = matchArgsByName(allDeps,analDeps,nameList) ;
3274 nameList.Delete() ;
3275 return result ;
3276}
3277
3278
3279
3280////////////////////////////////////////////////////////////////////////////////
3281/// Utility function for use in getAnalyticalIntegral(). If the
3282/// contents of proxies a,b,c,d occur in set 'allDeps' then the arguments
3283/// held in a,b,c,d are copied from allDeps to analDeps
3284
3285bool RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
3286 const RooArgProxy& a, const RooArgProxy& b,
3287 const RooArgProxy& c, const RooArgProxy& d) const
3288{
3289 TList nameList ;
3290 nameList.Add(new TObjString(a.absArg()->GetName())) ;
3291 nameList.Add(new TObjString(b.absArg()->GetName())) ;
3292 nameList.Add(new TObjString(c.absArg()->GetName())) ;
3293 nameList.Add(new TObjString(d.absArg()->GetName())) ;
3294 bool result = matchArgsByName(allDeps,analDeps,nameList) ;
3295 nameList.Delete() ;
3296 return result ;
3297}
3298
3299
3300////////////////////////////////////////////////////////////////////////////////
3301/// Utility function for use in getAnalyticalIntegral(). If the
3302/// contents of 'refset' occur in set 'allDeps' then the arguments
3303/// held in 'refset' are copied from allDeps to analDeps.
3304
3305bool RooAbsReal::matchArgs(const RooArgSet& allDeps, RooArgSet& analDeps,
3306 const RooArgSet& refset) const
3307{
3308 TList nameList ;
3309 for(RooAbsArg * arg : refset) {
3310 nameList.Add(new TObjString(arg->GetName())) ;
3311 }
3312
3313 bool result = matchArgsByName(allDeps,analDeps,nameList) ;
3314 nameList.Delete() ;
3315 return result ;
3316}
3317
3318
3319
3320////////////////////////////////////////////////////////////////////////////////
3321/// Check if allArgs contains matching elements for each name in nameList. If it does,
3322/// add the corresponding args from allArgs to matchedArgs and return true. Otherwise
3323/// return false and do not change matchedArgs.
3324
3325bool RooAbsReal::matchArgsByName(const RooArgSet &allArgs, RooArgSet &matchedArgs,
3326 const TList &nameList) const
3327{
3328 RooArgSet matched("matched");
3329 bool isMatched(true);
3330 for(auto * name : static_range_cast<TObjString*>(nameList)) {
3331 RooAbsArg *found= allArgs.find(name->String().Data());
3332 if(found) {
3333 matched.add(*found);
3334 }
3335 else {
3336 isMatched= false;
3337 break;
3338 }
3339 }
3340
3341 // nameList may not contain multiple entries with the same name
3342 // that are both matched
3343 if (isMatched && int(matched.size())!=nameList.GetSize()) {
3344 isMatched = false ;
3345 }
3346
3347 if(isMatched) matchedArgs.add(matched);
3348 return isMatched;
3349}
3350
3351
3352
3353////////////////////////////////////////////////////////////////////////////////
3354/// Returns the default numeric integration configuration for all RooAbsReals
3355
3357{
3359}
3360
3361
3362////////////////////////////////////////////////////////////////////////////////
3363/// Returns the specialized integrator configuration for _this_ RooAbsReal.
3364/// If this object has no specialized configuration, a null pointer is returned.
3365
3367{
3368 return _specIntegratorConfig.get();
3369}
3370
3371
3372////////////////////////////////////////////////////////////////////////////////
3373/// Returns the specialized integrator configuration for _this_ RooAbsReal.
3374/// If this object has no specialized configuration, a null pointer is returned,
3375/// unless createOnTheFly is true in which case a clone of the default integrator
3376/// configuration is created, installed as specialized configuration, and returned
3377
3379{
3380 if (!_specIntegratorConfig && createOnTheFly) {
3381 _specIntegratorConfig = std::make_unique<RooNumIntConfig>(*defaultIntegratorConfig()) ;
3382 }
3383 return _specIntegratorConfig.get();
3384}
3385
3386
3387
3388////////////////////////////////////////////////////////////////////////////////
3389/// Return the numeric integration configuration used for this object. If
3390/// a specialized configuration was associated with this object, that configuration
3391/// is returned, otherwise the default configuration for all RooAbsReals is returned
3392
3394{
3395 const RooNumIntConfig* config = specialIntegratorConfig() ;
3396 if (config) return config ;
3397 return defaultIntegratorConfig() ;
3398}
3399
3400
3401////////////////////////////////////////////////////////////////////////////////
3402/// Return the numeric integration configuration used for this object. If
3403/// a specialized configuration was associated with this object, that configuration
3404/// is returned, otherwise the default configuration for all RooAbsReals is returned
3405
3407{
3409 if (config) return config ;
3410 return defaultIntegratorConfig() ;
3411}
3412
3413
3414
3415////////////////////////////////////////////////////////////////////////////////
3416/// Set the given integrator configuration as default numeric integration
3417/// configuration for this object
3418
3420{
3421 _specIntegratorConfig = std::make_unique<RooNumIntConfig>(config);
3422}
3423
3424
3425
3426////////////////////////////////////////////////////////////////////////////////
3427/// Remove the specialized numeric integration configuration associated
3428/// with this object
3429
3431{
3432 _specIntegratorConfig.reset();
3433}
3434
3435
3436
3437
3438////////////////////////////////////////////////////////////////////////////////
3439/// Interface function to force use of a given set of observables
3440/// to interpret function value. Needed for functions or p.d.f.s
3441/// whose shape depends on the choice of normalization such as
3442/// RooAddPdf
3443
3445{
3446}
3447
3448
3449
3450
3451////////////////////////////////////////////////////////////////////////////////
3452/// Interface function to force use of a given normalization range
3453/// to interpret function value. Needed for functions or p.d.f.s
3454/// whose shape depends on the choice of normalization such as
3455/// RooAddPdf
3456
3458{
3459}
3460
3461
3462
3463////////////////////////////////////////////////////////////////////////////////
3464/// Advertise capability to determine maximum value of function for given set of
3465/// observables. If no direct generator method is provided, this information
3466/// will assist the accept/reject generator to operate more efficiently as
3467/// it can skip the initial trial sampling phase to empirically find the function
3468/// maximum
3469
3471{
3472 return 0 ;
3473}
3474
3475
3476
3477////////////////////////////////////////////////////////////////////////////////
3478/// Return maximum value for set of observables identified by code assigned
3479/// in getMaxVal
3480
3481double RooAbsReal::maxVal(Int_t /*code*/) const
3482{
3483 assert(1) ;
3484 return 0 ;
3485}
3486
3487
3488
3489////////////////////////////////////////////////////////////////////////////////
3490/// Interface to insert remote error logging messages received by RooRealMPFE into current error logging stream.
3491
3492void RooAbsReal::logEvalError(const RooAbsReal* originator, const char* origName, const char* message, const char* serverValueString)
3493{
3494 if (_evalErrorMode==Ignore) {
3495 return ;
3496 }
3497
3499 _evalErrorCount++ ;
3500 return ;
3501 }
3502
3503 static bool inLogEvalError = false ;
3504
3505 if (inLogEvalError) {
3506 return ;
3507 }
3508 inLogEvalError = true ;
3509
3510 EvalError ee ;
3511 ee.setMessage(message) ;
3512
3513 if (serverValueString) {
3514 ee.setServerValues(serverValueString) ;
3515 }
3516
3518 oocoutE(nullptr,Eval) << "RooAbsReal::logEvalError(" << "<STATIC>" << ") evaluation error, " << std::endl
3519 << " origin : " << origName << std::endl
3520 << " message : " << ee._msg << std::endl
3521 << " server values: " << ee._srvval << std::endl ;
3522 } else if (_evalErrorMode==CollectErrors) {
3523 _evalErrorList[originator].first = origName ;
3524 _evalErrorList[originator].second.push_back(ee) ;
3525 }
3526
3527
3528 inLogEvalError = false ;
3529}
3530
3531
3532
3533////////////////////////////////////////////////////////////////////////////////
3534/// Log evaluation error message. Evaluation errors may be routed through a different
3535/// protocol than generic RooFit warning message (which go straight through RooMsgService)
3536/// because evaluation errors can occur in very large numbers in the use of likelihood
3537/// evaluations. In logEvalError mode, controlled by global method enableEvalErrorLogging()
3538/// messages reported through this function are not printed but all stored in a list,
3539/// along with server values at the time of reporting. Error messages logged in this
3540/// way can be printed in a structured way, eliminating duplicates and with the ability
3541/// to truncate the list by printEvalErrors. This is the standard mode of error logging
3542/// during MINUIT operations. If enableEvalErrorLogging() is false, all errors
3543/// reported through this method are passed for immediate printing through RooMsgService.
3544/// A string with server names and values is constructed automatically for error logging
3545/// purposes, unless a custom string with similar information is passed as argument.
3546
3547void RooAbsReal::logEvalError(const char* message, const char* serverValueString) const
3548{
3549 if (_evalErrorMode==Ignore) {
3550 return ;
3551 }
3552
3554 _evalErrorCount++ ;
3555 return ;
3556 }
3557
3558 static bool inLogEvalError = false ;
3559
3560 if (inLogEvalError) {
3561 return ;
3562 }
3563 inLogEvalError = true ;
3564
3565 EvalError ee ;
3566 ee.setMessage(message) ;
3567
3568 if (serverValueString) {
3569 ee.setServerValues(serverValueString) ;
3570 } else {
3571 std::string srvval ;
3572 std::ostringstream oss ;
3573 bool first(true) ;
3574 for (Int_t i=0 ; i<numProxies() ; i++) {
3575 RooAbsProxy* p = getProxy(i) ;
3576 if (!p) continue ;
3577 //if (p->name()[0]=='!') continue ;
3578 if (first) {
3579 first=false ;
3580 } else {
3581 oss << ", " ;
3582 }
3583 p->print(oss,true) ;
3584 }
3585 ee.setServerValues(oss.str().c_str()) ;
3586 }
3587
3588 std::ostringstream oss2 ;
3590
3592 coutE(Eval) << "RooAbsReal::logEvalError(" << GetName() << ") evaluation error, " << std::endl
3593 << " origin : " << oss2.str() << std::endl
3594 << " message : " << ee._msg << std::endl
3595 << " server values: " << ee._srvval << std::endl ;
3596 } else if (_evalErrorMode==CollectErrors) {
3597 if (_evalErrorList[this].second.size() >= 2048) {
3598 // avoid overflowing the error list, so if there are very many, print
3599 // the oldest one first, and pop it off the list
3600 const EvalError& oee = _evalErrorList[this].second.front();
3601 // print to debug stream, since these would normally be suppressed, and
3602 // we do not want to increase the error count in the message service...
3603 ccoutD(Eval) << "RooAbsReal::logEvalError(" << GetName()
3604 << ") delayed evaluation error, " << std::endl
3605 << " origin : " << oss2.str() << std::endl
3606 << " message : " << oee._msg << std::endl
3607 << " server values: " << oee._srvval << std::endl ;
3608 _evalErrorList[this].second.pop_front();
3609 }
3610 _evalErrorList[this].first = oss2.str() ;
3611 _evalErrorList[this].second.push_back(ee) ;
3612 }
3613
3614 inLogEvalError = false ;
3615 //coutE(Tracing) << "RooAbsReal::logEvalError(" << GetName() << ") message = " << message << std::endl ;
3616}
3617
3618
3619
3620
3621////////////////////////////////////////////////////////////////////////////////
3622/// Clear the stack of evaluation error messages
3623
3625{
3627 return ;
3628 } else if (_evalErrorMode==CollectErrors) {
3629 _evalErrorList.clear() ;
3630 } else {
3631 _evalErrorCount = 0 ;
3632 }
3633}
3634
3635
3636////////////////////////////////////////////////////////////////////////////////
3637/// Retrieve bin boundaries if this distribution is binned in `obs`.
3638/// \param[in] obs Observable to retrieve boundaries for.
3639/// \param[in] xlo Beginning of range.
3640/// \param[in] xhi End of range.
3641/// \return The caller owns the returned std::list.
3642std::list<double>* RooAbsReal::binBoundaries(RooAbsRealLValue& /*obs*/, double /*xlo*/, double /*xhi*/) const {
3643 return nullptr;
3644}
3645
3646
3647////////////////////////////////////////////////////////////////////////////////
3648/// Interface for returning an optional hint for initial sampling points when constructing a curve projected on observable `obs`.
3649/// \param[in] obs Observable to retrieve sampling hint for.
3650/// \param[in] xlo Beginning of range.
3651/// \param[in] xhi End of range.
3652/// \return The caller owns the returned std::list.
3653std::list<double>* RooAbsReal::plotSamplingHint(RooAbsRealLValue& /*obs*/, double /*xlo*/, double /*xhi*/) const {
3654 return nullptr;
3655}
3656
3657////////////////////////////////////////////////////////////////////////////////
3658/// Print all outstanding logged evaluation error on the given ostream. If maxPerNode
3659/// is zero, only the number of errors for each source (object with unique name) is listed.
3660/// If maxPerNode is greater than zero, up to maxPerNode detailed error messages are shown
3661/// per source of errors. A truncation message is shown if there were more errors logged
3662/// than shown.
3663
3664void RooAbsReal::printEvalErrors(std::ostream& os, Int_t maxPerNode)
3665{
3666 if (_evalErrorMode == CountErrors) {
3667 os << _evalErrorCount << " errors counted" << std::endl ;
3668 }
3669
3670 if (maxPerNode<0) return ;
3671
3672 for(auto iter = _evalErrorList.begin();iter!=_evalErrorList.end() ; ++iter) {
3673 if (maxPerNode==0) {
3674
3675 // Only print node name with total number of errors
3676 os << iter->second.first ;
3677 //iter->first->printStream(os,kName|kClassName|kArgs,kInline) ;
3678 os << " has " << iter->second.second.size() << " errors" << std::endl ;
3679
3680 } else {
3681
3682 // Print node name and details of 'maxPerNode' errors
3683 os << iter->second.first << std::endl ;
3684 //iter->first->printStream(os,kName|kClassName|kArgs,kSingleLine) ;
3685
3686 Int_t i(0) ;
3687 for(auto iter2 = iter->second.second.begin();iter2!=iter->second.second.end() ; ++iter2, i++) {
3688 os << " " << iter2->_msg << " @ " << iter2->_srvval << std::endl ;
3689 if (i>maxPerNode) {
3690 os << " ... (remaining " << iter->second.second.size() - maxPerNode << " messages suppressed)" << std::endl ;
3691 break ;
3692 }
3693 }
3694 }
3695 }
3696}
3697
3698
3699
3700////////////////////////////////////////////////////////////////////////////////
3701/// Return the number of logged evaluation errors since the last clearing.
3702
3704{
3706 return _evalErrorCount ;
3707 }
3708
3709 Int_t ntot(0) ;
3710 for(auto const& elem : _evalErrorList) {
3711 ntot += elem.second.second.size() ;
3712 }
3713 return ntot ;
3714}
3715
3716
3717
3718////////////////////////////////////////////////////////////////////////////////
3719/// Fix the interpretation of the coefficient of any RooAddPdf component in
3720/// the expression tree headed by this object to the given set of observables.
3721///
3722/// If the force flag is false, the normalization choice is only fixed for those
3723/// RooAddPdf components that have the default 'automatic' interpretation of
3724/// coefficients (i.e. the interpretation is defined by the observables passed
3725/// to getVal()). If force is true, also RooAddPdf that already have a fixed
3726/// interpretation are changed to a new fixed interpretation.
3727
3728void RooAbsReal::fixAddCoefNormalization(const RooArgSet& addNormSet, bool force)
3729{
3730 std::unique_ptr<RooArgSet> compSet{getComponents()};
3731 for(auto * pdf : dynamic_range_cast<RooAbsPdf*>(*compSet)) {
3732 if (pdf) {
3733 if (!addNormSet.empty()) {
3734 pdf->selectNormalization(&addNormSet,force) ;
3735 } else {
3736 pdf->selectNormalization(nullptr,force) ;
3737 }
3738 }
3739 }
3740}
3741
3742
3743
3744////////////////////////////////////////////////////////////////////////////////
3745/// Fix the interpretation of the coefficient of any RooAddPdf component in
3746/// the expression tree headed by this object to the given set of observables.
3747///
3748/// If the force flag is false, the normalization range choice is only fixed for those
3749/// RooAddPdf components that currently use the default full domain to interpret their
3750/// coefficients. If force is true, also RooAddPdf that already have a fixed
3751/// interpretation range are changed to a new fixed interpretation range.
3752
3753void RooAbsReal::fixAddCoefRange(const char* rangeName, bool force)
3754{
3755 std::unique_ptr<RooArgSet> compSet{getComponents()};
3756 for(auto * pdf : dynamic_range_cast<RooAbsPdf*>(*compSet)) {
3757 if (pdf) {
3758 pdf->selectNormalizationRange(rangeName,force) ;
3759 }
3760 }
3761}
3762
3763
3764
3765////////////////////////////////////////////////////////////////////////////////
3766/// Interface method for function objects to indicate their preferred order of observables
3767/// for scanning their values into a (multi-dimensional) histogram or RooDataSet. The observables
3768/// to be ordered are offered in argument 'obs' and should be copied in their preferred
3769/// order into argument 'orderedObs', This default implementation indicates no preference
3770/// and copies the original order of 'obs' into 'orderedObs'
3771
3773{
3774 // Dummy implementation, do nothing
3775 orderedObs.removeAll() ;
3776 orderedObs.add(obs) ;
3777}
3778
3779
3780
3781////////////////////////////////////////////////////////////////////////////////
3782/// Calls createRunningIntegral(const RooArgSet&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&)
3783
3785{
3786 return createRunningIntegral(iset,RooFit::SupNormSet(nset)) ;
3787}
3788
3789
3790
3791////////////////////////////////////////////////////////////////////////////////
3792/// Create an object that represents the running integral of the function over one or more observables listed in iset, i.e.
3793/// \f[
3794/// \int_{x_\mathrm{lo}}^x f(x') \, \mathrm{d}x'
3795/// \f]
3796///
3797/// The actual integration calculation is only performed when the return object is evaluated. The name
3798/// of the integral object is automatically constructed from the name of the input function, the variables
3799/// it integrates and the range integrates over. The default strategy to calculate the running integrals is
3800///
3801/// - If the integrand (this object) supports analytical integration, construct an integral object
3802/// that calculate the running integrals value by calculating the analytical integral each
3803/// time the running integral object is evaluated
3804///
3805/// - If the integrand (this object) requires numeric integration to construct the running integral
3806/// create an object of class RooNumRunningInt which first samples the entire function and integrates
3807/// the sampled function numerically. This method has superior performance as there is no need to
3808/// perform a full (numeric) integration for each evaluation of the running integral object, but
3809/// only when one of its parameters has changed.
3810///
3811/// The choice of strategy can be changed with the ScanAll() argument, which forces the use of the
3812/// scanning technique implemented in RooNumRunningInt for all use cases, and with the ScanNone()
3813/// argument which forces the 'integrate each evaluation' technique for all use cases. The sampling
3814/// granularity for the scanning technique can be controlled with the ScanParameters technique
3815/// which allows to specify the number of samples to be taken, and to which order the resulting
3816/// running integral should be interpolated. The default values are 1000 samples and 2nd order
3817/// interpolation.
3818///
3819/// The following named arguments are accepted
3820/// | | Effect on integral creation
3821/// |-|-------------------------------
3822/// | `SupNormSet(const RooArgSet&)` | Observables over which should be normalized _in addition_ to the integration observables
3823/// | `ScanParameters(Int_t nbins, Int_t intOrder)` | Parameters for scanning technique of making CDF: number of sampled bins and order of interpolation applied on numeric cdf
3824/// | `ScanNum()` | Apply scanning technique if cdf integral involves numeric integration
3825/// | `ScanAll()` | Always apply scanning technique
3826/// | `ScanNone()` | Never apply scanning technique
3827
3829 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
3830 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
3831{
3832 // Define configuration for this method
3833 RooCmdConfig pc("RooAbsReal::createRunningIntegral(" + std::string(GetName()) + ")");
3834 pc.defineSet("supNormSet","SupNormSet",0,nullptr) ;
3835 pc.defineInt("numScanBins","ScanParameters",0,1000) ;
3836 pc.defineInt("intOrder","ScanParameters",1,2) ;
3837 pc.defineInt("doScanNum","ScanNum",0,1) ;
3838 pc.defineInt("doScanAll","ScanAll",0,0) ;
3839 pc.defineInt("doScanNon","ScanNone",0,0) ;
3840 pc.defineMutex("ScanNum","ScanAll","ScanNone") ;
3841
3842 // Process & check varargs
3843 pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
3844 if (!pc.ok(true)) {
3845 return nullptr ;
3846 }
3847
3848 // Extract values from named arguments
3849 RooArgSet nset ;
3850 if (const RooArgSet* snset = pc.getSet("supNormSet",nullptr)) {
3851 nset.add(*snset) ;
3852 }
3853 Int_t numScanBins = pc.getInt("numScanBins") ;
3854 Int_t intOrder = pc.getInt("intOrder") ;
3855 Int_t doScanNum = pc.getInt("doScanNum") ;
3856 Int_t doScanAll = pc.getInt("doScanAll") ;
3857 Int_t doScanNon = pc.getInt("doScanNon") ;
3858
3859 // If scanning technique is not requested make integral-based cdf and return
3860 if (doScanNon) {
3861 return createIntRI(iset,nset) ;
3862 }
3863 if (doScanAll) {
3864 return createScanRI(iset,nset,numScanBins,intOrder) ;
3865 }
3866 if (doScanNum) {
3867 std::unique_ptr<RooAbsReal> tmp{createIntegral(iset)} ;
3868 Int_t isNum= !static_cast<RooRealIntegral&>(*tmp).numIntRealVars().empty();
3869
3870 if (isNum) {
3871 coutI(NumIntegration) << "RooAbsPdf::createRunningIntegral(" << GetName() << ") integration over observable(s) " << iset << " involves numeric integration," << std::endl
3872 << " constructing cdf though numeric integration of sampled pdf in " << numScanBins << " bins and applying order "
3873 << intOrder << " interpolation on integrated histogram." << std::endl
3874 << " To override this choice of technique use argument ScanNone(), to change scan parameters use ScanParameters(nbins,order) argument" << std::endl ;
3875 }
3876
3877 return isNum ? createScanRI(iset,nset,numScanBins,intOrder) : createIntRI(iset,nset) ;
3878 }
3879 return nullptr;
3880}
3881
3882
3883
3884////////////////////////////////////////////////////////////////////////////////
3885/// Utility function for createRunningIntegral that construct an object
3886/// implementing the numeric scanning technique for calculating the running integral
3887
3889{
3890 std::string name = std::string(GetName()) + "_NUMRUNINT_" + integralNameSuffix(iset,&nset).Data() ;
3891 RooRealVar* ivar = static_cast<RooRealVar*>(iset.first()) ;
3892 ivar->setBins(numScanBins,"numcdf") ;
3893 auto ret = std::make_unique<RooNumRunningInt>(name.c_str(),name.c_str(),*this,*ivar,"numrunint") ;
3894 ret->setInterpolationOrder(intOrder) ;
3895 return RooFit::makeOwningPtr<RooAbsReal>(std::move(ret));
3896}
3897
3898
3899
3900////////////////////////////////////////////////////////////////////////////////
3901/// Utility function for createRunningIntegral. It creates an
3902/// object implementing the standard (analytical) integration
3903/// technique for calculating the running integral.
3904
3906{
3907 // Make list of input arguments keeping only RooRealVars
3908 RooArgList ilist ;
3909 for(RooAbsArg * arg : iset) {
3910 if (dynamic_cast<RooRealVar*>(arg)) {
3911 ilist.add(*arg) ;
3912 } else {
3913 coutW(InputArguments) << "RooAbsPdf::createRunningIntegral(" << GetName() << ") WARNING ignoring non-RooRealVar input argument " << arg->GetName() << std::endl ;
3914 }
3915 }
3916
3917 RooArgList cloneList ;
3918 RooArgList loList ;
3919 RooArgSet clonedBranchNodes ;
3920
3921 // Setup customizer that stores all cloned branches in our non-owning list
3922 RooCustomizer cust(*this,"cdf") ;
3923 cust.setCloneBranchSet(clonedBranchNodes) ;
3924 cust.setOwning(false) ;
3925
3926 // Make integration observable x_prime for each observable x as well as an x_lowbound
3927 for(auto * rrv : static_range_cast<RooRealVar*>(ilist)) {
3928
3929 // Make clone x_prime of each c.d.f observable x represening running integral
3930 RooRealVar* cloneArg = static_cast<RooRealVar*>(rrv->clone(Form("%s_prime",rrv->GetName()))) ;
3931 cloneList.add(*cloneArg) ;
3932 cust.replaceArg(*rrv,*cloneArg) ;
3933
3934 // Make clone x_lowbound of each c.d.f observable representing low bound of x
3935 RooRealVar* cloneLo = static_cast<RooRealVar*>(rrv->clone(Form("%s_lowbound",rrv->GetName()))) ;
3936 cloneLo->setVal(rrv->getMin()) ;
3937 loList.add(*cloneLo) ;
3938
3939 // Make parameterized binning from [x_lowbound,x] for each x_prime
3940 RooParamBinning pb(*cloneLo,*rrv,100) ;
3941 cloneArg->setBinning(pb,"CDF") ;
3942
3943 }
3944
3945 RooAbsReal* tmp = static_cast<RooAbsReal*>(cust.build()) ;
3946
3947 // Construct final normalization set for c.d.f = integrated observables + any extra specified by user
3948 RooArgSet finalNset(nset) ;
3949 finalNset.add(cloneList,true) ;
3950 std::unique_ptr<RooAbsReal> cdf{tmp->createIntegral(cloneList,finalNset,"CDF")};
3951
3952 // Transfer ownership of cloned items to top-level c.d.f object
3953 cdf->addOwnedComponents(*tmp) ;
3954 cdf->addOwnedComponents(cloneList) ;
3955 cdf->addOwnedComponents(loList) ;
3956
3957 return RooFit::makeOwningPtr(std::move(cdf));
3958}
3959
3960
3961////////////////////////////////////////////////////////////////////////////////
3962/// Return a RooFunctor object bound to this RooAbsReal with given definition of observables
3963/// and parameters
3964
3965RooFunctor* RooAbsReal::functor(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
3966{
3967 RooArgSet realObs;
3968 getObservables(&obs, realObs);
3969 if (realObs.size() != obs.size()) {
3970 coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified observables are not variables of this p.d.f" << std::endl ;
3971 return nullptr;
3972 }
3973 RooArgSet realPars;
3974 getObservables(&pars, realPars);
3975 if (realPars.size() != pars.size()) {
3976 coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified parameters are not variables of this p.d.f" << std::endl ;
3977 return nullptr;
3978 }
3979
3980 return new RooFunctor(*this,obs,pars,nset) ;
3981}
3982
3983
3984
3985////////////////////////////////////////////////////////////////////////////////
3986/// Return a ROOT TF1,2,3 object bound to this RooAbsReal with given definition of observables
3987/// and parameters
3988
3989TF1* RooAbsReal::asTF(const RooArgList& obs, const RooArgList& pars, const RooArgSet& nset) const
3990{
3991 // Check that specified input are indeed variables of this function
3992 RooArgSet realObs;
3993 getObservables(&obs, realObs) ;
3994 if (realObs.size() != obs.size()) {
3995 coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified observables are not variables of this p.d.f" << std::endl ;
3996 return nullptr ;
3997 }
3998 RooArgSet realPars;
3999 getObservables(&pars, realPars) ;
4000 if (realPars.size() != pars.size()) {
4001 coutE(InputArguments) << "RooAbsReal::functor(" << GetName() << ") ERROR: one or more specified parameters are not variables of this p.d.f" << std::endl ;
4002 return nullptr ;
4003 }
4004
4005 // Check that all obs and par are of type RooRealVar
4006 for (std::size_t i=0 ; i<obs.size() ; i++) {
4007 if (dynamic_cast<RooRealVar*>(obs.at(i))==nullptr) {
4008 coutE(ObjectHandling) << "RooAbsReal::asTF(" << GetName() << ") ERROR: proposed observable " << obs.at(0)->GetName() << " is not of type RooRealVar" << std::endl ;
4009 return nullptr ;
4010 }
4011 }
4012 for (std::size_t i=0 ; i<pars.size() ; i++) {
4013 if (dynamic_cast<RooRealVar*>(pars.at(i))==nullptr) {
4014 coutE(ObjectHandling) << "RooAbsReal::asTF(" << GetName() << ") ERROR: proposed parameter " << pars.at(0)->GetName() << " is not of type RooRealVar" << std::endl ;
4015 return nullptr ;
4016 }
4017 }
4018
4019 // Create functor and TFx of matching dimension
4020 TF1* tf=nullptr ;
4021 RooFunctor* f ;
4022 switch(obs.size()) {
4023 case 1: {
4024 RooRealVar* x = static_cast<RooRealVar*>(obs.at(0)) ;
4025 f = functor(obs,pars,nset) ;
4026 tf = new TF1(GetName(),f,x->getMin(),x->getMax(),pars.size()) ;
4027 break ;
4028 }
4029 case 2: {
4030 RooRealVar* x = static_cast<RooRealVar*>(obs.at(0)) ;
4031 RooRealVar* y = static_cast<RooRealVar*>(obs.at(1)) ;
4032 f = functor(obs,pars,nset) ;
4033 tf = new TF2(GetName(),f,x->getMin(),x->getMax(),y->getMin(),y->getMax(),pars.size()) ;
4034 break ;
4035 }
4036 case 3: {
4037 RooRealVar* x = static_cast<RooRealVar*>(obs.at(0)) ;
4038 RooRealVar* y = static_cast<RooRealVar*>(obs.at(1)) ;
4039 RooRealVar* z = static_cast<RooRealVar*>(obs.at(2)) ;
4040 f = functor(obs,pars,nset) ;
4041 tf = new TF3(GetName(),f,x->getMin(),x->getMax(),y->getMin(),y->getMax(),z->getMin(),z->getMax(),pars.size()) ;
4042 break ;
4043 }
4044 default:
4045 coutE(InputArguments) << "RooAbsReal::asTF(" << GetName() << ") ERROR: " << obs.size()
4046 << " observables specified, but a ROOT TFx can only have 1,2 or 3 observables" << std::endl ;
4047 return nullptr ;
4048 }
4049
4050 // Set initial parameter values of TFx to those of RooRealVars
4051 for (std::size_t i=0 ; i<pars.size() ; i++) {
4052 RooRealVar* p = static_cast<RooRealVar*>(pars.at(i)) ;
4053 tf->SetParameter(i,p->getVal()) ;
4054 tf->SetParName(i,p->GetName()) ;
4055 //tf->SetParLimits(i,p->getMin(),p->getMax()) ;
4056 }
4057
4058 return tf ;
4059}
4060
4061
4062////////////////////////////////////////////////////////////////////////////////
4063/// Return function representing first, second or third order derivative of this function
4064
4066{
4067 std::string name=Form("%s_DERIV_%s",GetName(),obs.GetName()) ;
4068 std::string title=Form("Derivative of %s w.r.t %s ",GetName(),obs.GetName()) ;
4069 return new RooDerivative(name.c_str(),title.c_str(),*this,obs,order,eps) ;
4070}
4071
4072
4073
4074////////////////////////////////////////////////////////////////////////////////
4075/// Return function representing first, second or third order derivative of this function
4076
4077RooDerivative* RooAbsReal::derivative(RooRealVar& obs, const RooArgSet& normSet, Int_t order, double eps)
4078{
4079 std::string name=Form("%s_DERIV_%s",GetName(),obs.GetName()) ;
4080 std::string title=Form("Derivative of %s w.r.t %s ",GetName(),obs.GetName()) ;
4081 return new RooDerivative(name.c_str(),title.c_str(),*this,obs,normSet,order,eps) ;
4082}
4083
4084
4085
4086////////////////////////////////////////////////////////////////////////////////
4087/// Return function representing moment of function of given order.
4088/// \param[in] obs Observable to calculate the moments for
4089/// \param[in] order Order of the moment
4090/// \param[in] central If true, the central moment is given by \f$ \langle (x- \langle x \rangle )^2 \rangle \f$
4091/// \param[in] takeRoot Calculate the square root
4092
4093RooAbsMoment* RooAbsReal::moment(RooRealVar& obs, Int_t order, bool central, bool takeRoot)
4094{
4095 std::string name=Form("%s_MOMENT_%d%s_%s",GetName(),order,(central?"C":""),obs.GetName()) ;
4096 std::string title=Form("%sMoment of order %d of %s w.r.t %s ",(central?"Central ":""),order,GetName(),obs.GetName()) ;
4097 if (order==1) return new RooFirstMoment(name.c_str(),title.c_str(),*this,obs) ;
4098 if (order==2) return new RooSecondMoment(name.c_str(),title.c_str(),*this,obs,central,takeRoot) ;
4099 return new RooMoment(name.c_str(),title.c_str(),*this,obs,order,central,takeRoot) ;
4100}
4101
4102
4103////////////////////////////////////////////////////////////////////////////////
4104/// Return function representing moment of p.d.f (normalized w.r.t given observables) of given order.
4105/// \param[in] obs Observable to calculate the moments for
4106/// \param[in] normObs Normalise w.r.t. these observables
4107/// \param[in] order Order of the moment
4108/// \param[in] central If true, the central moment is given by \f$ \langle (x- \langle x \rangle )^2 \rangle \f$
4109/// \param[in] takeRoot Calculate the square root
4110/// \param[in] intNormObs If true, the moment of the function integrated over all normalization observables is returned.
4111
4112RooAbsMoment* RooAbsReal::moment(RooRealVar& obs, const RooArgSet& normObs, Int_t order, bool central, bool takeRoot, bool intNormObs)
4113{
4114 std::string name=Form("%s_MOMENT_%d%s_%s",GetName(),order,(central?"C":""),obs.GetName()) ;
4115 std::string title=Form("%sMoment of order %d of %s w.r.t %s ",(central?"Central ":""),order,GetName(),obs.GetName()) ;
4116
4117 if (order==1) return new RooFirstMoment(name.c_str(),title.c_str(),*this,obs,normObs,intNormObs) ;
4118 if (order==2) return new RooSecondMoment(name.c_str(),title.c_str(),*this,obs,normObs,central,takeRoot,intNormObs) ;
4119 return new RooMoment(name.c_str(),title.c_str(),*this,obs,normObs,order,central,takeRoot,intNormObs) ;
4120}
4121
4122
4123
4124////////////////////////////////////////////////////////////////////////////////
4125///
4126/// Return value of x (in range xmin,xmax) at which function equals yval.
4127/// (Calculation is performed with Brent root finding algorithm)
4128
4129double RooAbsReal::findRoot(RooRealVar& x, double xmin, double xmax, double yval)
4130{
4131 double result(0) ;
4133 return result ;
4134}
4135
4136
4137
4138////////////////////////////////////////////////////////////////////////////////
4139/// Perform a \f$ \chi^2 \f$ fit to given histogram. By default the fit is executed through the MINUIT
4140/// commands MIGRAD, HESSE in succession
4141///
4142/// The following named arguments are supported
4143///
4144/// <table>
4145/// <tr><th> <th> Options to control construction of chi2
4146/// <tr><td> `Extended(bool flag)` <td> **Only applicable when fitting a RooAbsPdf**. Scale the normalized pdf by the number of events predicted by the model instead of scaling by the total data weight.
4147/// This imposes a constraint on the predicted number of events analogous to the extended term in a likelihood fit.
4148/// - If you don't pass this command, an extended fit will be done by default if the pdf makes a prediction on the number of events
4149/// (in RooFit jargon, "if the pdf can be extended").
4150/// - Passing `Extended(true)` when the the pdf makes no prediction on the expected number of events will result in error messages,
4151/// and the chi2 will fall back to the total data weight to scale the normalized pdf.
4152/// - There are cases where the fit **must** be done in extended mode. This happens for example when you have a RooAddPdf
4153/// where the coefficients represent component yields. If the fit is not extended, these coefficients will not be
4154/// well-defined, as the RooAddPdf always normalizes itself. If you pass `Extended(false)` in such a case, an error will be
4155/// printed and you'll most likely get garbage results.
4156/// <tr><td> `Range(const char* name)` <td> Fit only data inside range with given name
4157/// <tr><td> `Range(double lo, double hi)` <td> Fit only data inside given range. A range named "fit" is created on the fly on all observables.
4158/// Multiple comma separated range names can be specified.
4159/// <tr><td> `NumCPU(int num)` <td> Parallelize NLL calculation on num CPUs
4160/// <tr><td> `Optimize(bool flag)` <td> Activate constant term optimization (on by default)
4161/// <tr><td> `IntegrateBins()` <td> Integrate PDF within each bin. This sets the desired precision.
4162///
4163/// <tr><th> <th> Options to control flow of fit procedure
4164/// <tr><td> `InitialHesse(bool flag)` <td> Flag controls if HESSE before MIGRAD as well, off by default
4165/// <tr><td> `Hesse(bool flag)` <td> Flag controls if HESSE is run after MIGRAD, on by default
4166/// <tr><td> `Minos(bool flag)` <td> Flag controls if MINOS is run after HESSE, on by default
4167/// <tr><td> `Minos(const RooArgSet& set)` <td> Only run MINOS on given subset of arguments
4168/// <tr><td> `Save(bool flag)` <td> Flag controls if RooFitResult object is produced and returned, off by default
4169/// <tr><td> `Strategy(Int_t flag)` <td> Set Minuit strategy (0 through 2, default is 1)
4170///
4171/// <tr><th> <th> Options to control informational output
4172/// <tr><td> `Verbose(bool flag)` <td> Flag controls if verbose output is printed (NLL, parameter changes during fit
4173/// <tr><td> `Timer(bool flag)` <td> Time CPU and wall clock consumption of fit steps, off by default
4174/// <tr><td> `PrintLevel(Int_t level)` <td> Set Minuit print level (-1 through 3, default is 1). At -1 all RooFit informational
4175/// messages are suppressed as well
4176/// <tr><td> `Warnings(bool flag)` <td> Enable or disable MINUIT warnings (enabled by default)
4177/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per likelihood evaluation. A negative
4178/// value suppress output completely, a zero value will only print the error count per p.d.f component,
4179/// a positive value is will print details of each error up to numErr messages per p.d.f component.
4180/// </table>
4181///
4182
4184 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
4185 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
4186{
4188 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
4189 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
4190 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
4191 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
4192 return chi2FitTo(data,l) ;
4193}
4194
4195
4196
4197////////////////////////////////////////////////////////////////////////////////
4198/// Calls RooAbsReal::createChi2(RooDataSet& data, const RooLinkedList& cmdList) and returns fit result.
4199///
4200/// List of possible commands in the `cmdList`:
4201///
4202/// <table>
4203/// <tr><th> Type of CmdArg <th> Effect on \f$ \chi^2 \f$
4204/// <tr><td>
4205/// <tr><td> `DataError()` <td> Choose between:
4206/// - RooAbsData::Expected: Expected Poisson error (\f$ \sqrt{n_\text{expected}} \f$ from the PDF).
4207/// - RooAbsData::SumW2: The observed error from the square root of the sum of weights,
4208/// i.e., symmetric errors calculated with the standard deviation of a Poisson distribution.
4209/// - RooAbsData::Poisson: Asymmetric errors from the central 68 % interval around a Poisson distribution with mean \f$ n_\text{observed} \f$.
4210/// If for a given bin \f$ n_\text{expected} \f$ is lower than the \f$ n_\text{observed} \f$, the lower uncertainty is taken
4211/// (e.g., the difference between the mean and the 16 % quantile).
4212/// If \f$ n_\text{expected} \f$ is higher than \f$ n_\text{observed} \f$, the higher uncertainty is taken
4213/// (e.g., the difference between the 84 % quantile and the mean).
4214/// - RooAbsData::Auto (default): RooAbsData::Expected for unweighted data, RooAbsData::SumW2 for weighted data.
4215/// <tr><td>
4216/// `Extended()` <td> Use expected number of events of an extended p.d.f as normalization
4217/// <tr><td>
4218/// NumCPU() <td> Activate parallel processing feature
4219/// <tr><td>
4220/// Range() <td> Calculate \f$ \chi^2 \f$ only in selected region
4221/// <tr><td>
4222/// Verbose() <td> Verbose output of GOF framework
4223/// <tr><td>
4224/// IntegrateBins() <td> Integrate PDF within each bin. This sets the desired precision. Only useful for binned fits.
4225/// <tr><td> `SumCoefRange()` <td> Set the range in which to interpret the coefficients of RooAddPdf components
4226/// <tr><td> `SplitRange()` <td> Fit ranges used in different categories get named after the category.
4227/// Using `Range("range"), SplitRange()` as switches, different ranges could be set like this:
4228/// ```
4229/// myVariable.setRange("range_pi0", 135, 210);
4230/// myVariable.setRange("range_gamma", 50, 210);
4231/// ```
4232/// <tr><td> `ConditionalObservables(Args_t &&... argsOrArgSet)` <td> Define projected observables.
4233/// Arguments can either be multiple RooRealVar or a single RooArgSet containing them.
4234///
4235/// </table>
4236
4238{
4239 return RooFit::makeOwningPtr(RooFit::FitHelpers::fitTo(*this, data, cmdList, true));
4240}
4241
4242
4243
4244
4245////////////////////////////////////////////////////////////////////////////////
4246/// Create a \f$ \chi^2 \f$ variable from a histogram and this function.
4247///
4248/// \param arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 ordered arguments
4249///
4250/// The list of supported command arguments is given in the documentation for
4251/// RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataHist& hdata, const RooCmdArg&,const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,const RooCmdArg&,const RooCmdArg&).
4252///
4253/// \param data Histogram with data
4254/// \return \f$ \chi^2 \f$ variable
4255
4257 const RooCmdArg &arg3, const RooCmdArg &arg4,
4258 const RooCmdArg &arg5, const RooCmdArg &arg6,
4259 const RooCmdArg &arg7, const RooCmdArg &arg8)
4260{
4262 l.Add((TObject *)&arg1);
4263 l.Add((TObject *)&arg2);
4264 l.Add((TObject *)&arg3);
4265 l.Add((TObject *)&arg4);
4266 l.Add((TObject *)&arg5);
4267 l.Add((TObject *)&arg6);
4268 l.Add((TObject *)&arg7);
4269 l.Add((TObject *)&arg8);
4270 return createChi2(data, l);
4271}
4272
4273////////////////////////////////////////////////////////////////////////////////
4274/// \see RooAbsReal::createChi2(RooDataHist&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4275/// \param data hist data
4276/// \param cmdList List with RooCmdArg() from the table
4277
4279{
4280 return RooFit::makeOwningPtr(RooFit::FitHelpers::createChi2(*this, data, cmdList));
4281}
4282
4283////////////////////////////////////////////////////////////////////////////////
4284/// Perform a 2-D \f$ \chi^2 \f$ fit using a series of x and y values stored in the dataset `xydata`.
4285/// The y values can either be the event weights, or can be another column designated
4286/// by the YVar() argument. The y value must have errors defined for the \f$ \chi^2 \f$ to
4287/// be well defined.
4288///
4289/// <table>
4290/// <tr><th><th> Options to control construction of the chi-square
4291/// <tr><td> `YVar(RooRealVar& yvar)` <td> Designate given column in dataset as Y value
4292/// <tr><td> `Integrate(bool flag)` <td> Integrate function over range specified by X errors
4293/// rather than take value at bin center.
4294///
4295/// <tr><th><th> Options to control flow of fit procedure
4296/// <tr><td> `InitialHesse(bool flag)` <td> Flag controls if HESSE before MIGRAD as well, off by default
4297/// <tr><td> `Hesse(bool flag)` <td> Flag controls if HESSE is run after MIGRAD, on by default
4298/// <tr><td> `Minos(bool flag)` <td> Flag controls if MINOS is run after HESSE, on by default
4299/// <tr><td> `Minos(const RooArgSet& set)` <td> Only run MINOS on given subset of arguments
4300/// <tr><td> `Save(bool flag)` <td> Flag controls if RooFitResult object is produced and returned, off by default
4301/// <tr><td> `Strategy(Int_t flag)` <td> Set Minuit strategy (0 through 2, default is 1)
4302///
4303/// <tr><th><th> Options to control informational output
4304/// <tr><td> `Verbose(bool flag)` <td> Flag controls if verbose output is printed (NLL, parameter changes during fit
4305/// <tr><td> `Timer(bool flag)` <td> Time CPU and wall clock consumption of fit steps, off by default
4306/// <tr><td> `PrintLevel(Int_t level)` <td> Set Minuit print level (-1 through 3, default is 1). At -1 all RooFit informational
4307/// messages are suppressed as well
4308/// <tr><td> `Warnings(bool flag)` <td> Enable or disable MINUIT warnings (enabled by default)
4309/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per likelihood evaluation. A negative
4310/// value suppress output completely, a zero value will only print the error count per p.d.f component,
4311/// a positive value is will print details of each error up to numErr messages per p.d.f component.
4312/// </table>
4313
4315 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
4316 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
4317{
4319 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
4320 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
4321 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
4322 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
4323 return chi2FitTo(xydata,l) ;
4324}
4325
4326
4327
4328
4329////////////////////////////////////////////////////////////////////////////////
4330/// \copydoc RooAbsReal::chi2FitTo(RooDataSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4331
4333{
4334 return RooFit::makeOwningPtr(RooFit::FitHelpers::fitTo(*this, xydata, cmdList, true));
4335}
4336
4337
4338
4339
4340////////////////////////////////////////////////////////////////////////////////
4341/// Create a \f$ \chi^2 \f$ from a series of x and y values stored in a dataset.
4342/// The y values can either be the event weights (default), or can be another column designated
4343/// by the YVar() argument. The y value must have errors defined for the \f$ \chi^2 \f$ to
4344/// be well defined.
4345///
4346/// The following named arguments are supported
4347///
4348/// | | Options to control construction of the \f$ \chi^2 \f$
4349/// |-|-----------------------------------------
4350/// | `YVar(RooRealVar& yvar)` | Designate given column in dataset as Y value
4351/// | `Integrate(bool flag)` | Integrate function over range specified by X errors rather than take value at bin center.
4352///
4353
4355 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
4356 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
4357{
4359 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
4360 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
4361 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
4362 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
4363 return createChi2(data,l) ;
4364}
4365
4366
4367////////////////////////////////////////////////////////////////////////////////
4368/// See RooAbsReal::createChi2(RooDataSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4369
4371{
4372 return RooFit::makeOwningPtr(RooFit::FitHelpers::createChi2(*this, data, cmdList));
4373}
4374
4375
4376
4377////////////////////////////////////////////////////////////////////////////////
4378/// Return current evaluation error logging mode.
4379
4381{
4382 return _evalErrorMode ;
4383}
4384
4385////////////////////////////////////////////////////////////////////////////////
4386/// Set evaluation error logging mode. Options are
4387///
4388/// PrintErrors - Print each error through RooMsgService() as it occurs
4389/// CollectErrors - Accumulate errors, but do not print them. A subsequent call
4390/// to printEvalErrors() will print a summary
4391/// CountErrors - Accumulate error count, but do not print them.
4392///
4393
4395{
4396 _evalErrorMode = m;
4397}
4398
4399
4400////////////////////////////////////////////////////////////////////////////////
4401
4403{
4404 std::string plist ;
4405 for (auto const* arg : paramVars) {
4406 if (!dependsOnValue(*arg)) {
4407 coutW(InputArguments) << "RooAbsReal::setParameterizeIntegral(" << GetName()
4408 << ") function does not depend on listed parameter " << arg->GetName() << ", ignoring" << std::endl ;
4409 continue ;
4410 }
4411 if (!plist.empty()) plist += ":" ;
4412 plist += arg->GetName() ;
4413 }
4414 setStringAttribute("CACHEPARAMINT",plist.c_str()) ;
4415}
4416
4417
4418/** Base function for computing multiple values of a RooAbsReal
4419\param output The array where the results are stored
4420\param nEvents The number of events to be processed
4421\param dataMap A std::map containing the input data for the computations
4422**/
4423void RooAbsReal::computeBatch(double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const {
4424 // Find all servers that are serving real numbers to us, retrieve their batch data,
4425 // and switch them into "always clean" operating mode, so they return always the last-set value.
4426 struct ServerData {
4427 RooAbsArg* server;
4428 std::span<const double> batch;
4429 double oldValue;
4430 RooAbsArg::OperMode oldOperMode;
4431 bool oldValueDirty;
4432 bool oldShapeDirty;
4433 };
4434 std::vector<ServerData> ourServers;
4435 ourServers.reserve(servers().size());
4436
4437 for (auto server : servers()) {
4438 auto serverValues = dataMap.at(server);
4439 if(serverValues.empty()) continue;
4440
4441 // maybe we are still missing inhibit dirty here
4442 auto oldOperMode = server->operMode();
4443 // See note at the bottom of this function to learn why we can only set
4444 // the operation mode to "always clean" if there are no other value
4445 // clients.
4446 server->setOperMode(RooAbsArg::AClean);
4447 ourServers.push_back({server,
4448 serverValues,
4449 server->isCategory() ? static_cast<RooAbsCategory const*>(server)->getCurrentIndex() : static_cast<RooAbsReal const*>(server)->_value,
4450 oldOperMode,
4451 server->_valueDirty,
4452 server->_shapeDirty});
4453 // Prevent the server from evaluating; just return cached result, which we will side load:
4454 }
4455
4456
4457 // Make sure that we restore all state when we finish:
4458 struct RestoreStateRAII {
4459 RestoreStateRAII(std::vector<ServerData>& servers) :
4460 _servers{servers} { }
4461
4462 ~RestoreStateRAII() {
4463 for (auto& serverData : _servers) {
4464 serverData.server->setCachedValue(serverData.oldValue, true);
4465 serverData.server->setOperMode(serverData.oldOperMode);
4466 serverData.server->_valueDirty = serverData.oldValueDirty;
4467 serverData.server->_shapeDirty = serverData.oldShapeDirty;
4468 }
4469 }
4470
4471 std::vector<ServerData>& _servers;
4472 } restoreState{ourServers};
4473
4474
4475 // Advising to implement the batch interface makes only sense if the batch was not a scalar.
4476 // Otherwise, there would be no speedup benefit.
4478 coutI(FastEvaluations) << "The class " << ClassName() << " does not implement the faster batch evaluation interface."
4479 << " Consider requesting or implementing it to benefit from a speed up." << std::endl;
4480 }
4481
4482
4483 // For each event, write temporary values into our servers' caches, and run a single-value computation.
4484
4485 for (std::size_t i=0; i < nEvents; ++i) {
4486 for (auto& serv : ourServers) {
4487 serv.server->setCachedValue(serv.batch[std::min(i, serv.batch.size()-1)], false);
4488 }
4489
4490 output[i] = evaluate();
4491 }
4492}
4493
4494////////////////////////////////////////////////////////////////////////////////
4495/// This function defines the analytical integral translation for the class.
4496///
4497/// \param[in] code The code that decides the integrands.
4498/// \param[in] rangeName Name of the normalization range.
4499/// \param[in] ctx An object to manage auxiliary information for code-squashing.
4500///
4501/// \returns The representative code string of the integral for the given object.
4502std::string RooAbsReal::buildCallToAnalyticIntegral(Int_t /* code */, const char * /* rangeName */,
4503 RooFit::Detail::CodeSquashContext & /*ctx*/) const
4504{
4505 std::stringstream errorMsg;
4506 errorMsg << "An analytical integral function for class \"" << ClassName() << "\" has not yet been implemented.";
4507 coutE(Minimization) << errorMsg.str() << std::endl;
4508 throw std::runtime_error(errorMsg.str().c_str());
4509}
4510
4511double RooAbsReal::_DEBUG_getVal(const RooArgSet* normalisationSet) const {
4512
4513 const bool tmpFast = _fast;
4514 const double tmp = _value;
4515
4516 double fullEval = 0.;
4517 try {
4518 fullEval = getValV(normalisationSet);
4519 }
4520 catch (CachingError& error) {
4521 throw CachingError(std::move(error),
4522 FormatPdfTree() << *this);
4523 }
4524
4525 const double ret = (_fast && !_inhibitDirty) ? _value : fullEval;
4526
4527 if (std::isfinite(ret) && ( ret != 0. ? (ret - fullEval)/ret : ret - fullEval) > 1.E-9) {
4528#ifndef NDEBUG
4530#endif
4531 FormatPdfTree formatter;
4532 formatter << "--> (Scalar computation wrong here:)\n"
4533 << GetName() << " " << this << " _fast=" << tmpFast
4534 << "\n\tcached _value=" << std::setprecision(16) << tmp
4535 << "\n\treturning =" << ret
4536 << "\n\trecomputed =" << fullEval
4537 << "\n\tnew _value =" << _value << "] ";
4538 formatter << "\nServers:";
4539 for (const auto server : _serverList) {
4540 formatter << "\n ";
4541 server->printStream(formatter.stream(), kName | kClassName | kArgs | kExtras | kAddress | kValue, kInline);
4542 }
4543
4544 throw CachingError(formatter);
4545 }
4546
4547 return ret;
4548}
4549
4550
4551bool RooAbsReal::redirectServersHook(const RooAbsCollection & newServerList, bool mustReplaceAll,
4552 bool nameChange, bool isRecursiveStep)
4553{
4555 return RooAbsArg::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursiveStep);
4556}
4557
4558
4559////////////////////////////////////////////////////////////////////////////////
4560
4562{
4563 for (RooAbsArg* arg : servers()) {
4564 if(auto realArg = dynamic_cast<RooAbsReal*>(arg)) {
4565 realArg->enableOffsetting(flag) ;
4566 }
4567 }
4568}
4569
4570
4571RooAbsReal::Ref::Ref(double val) : _ref{RooFit::RooConst(val)} {}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define cxcoutD(a)
#define coutW(a)
#define ccoutP(a)
#define dologD(a)
#define coutF(a)
#define oocoutE(o, a)
#define coutE(a)
#define ccoutW(a)
#define ccoutD(a)
int Int_t
Definition RtypesCore.h:45
float Size_t
Definition RtypesCore.h:96
char Text_t
Definition RtypesCore.h:62
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
static void indent(ostringstream &buf, int indent_level)
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t SetLineColor
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format
char name[80]
Definition TGX11.cxx:110
float xmin
float ymin
float xmax
float ymax
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2467
R__EXTERN TSystem * gSystem
Definition TSystem.h:560
std::ostream & stream()
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
TIterator Use servers() and begin()
bool dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr, bool valueOnly=false) const
Test whether we depend on (ie, are served by) any object in the specified collection.
bool recursiveRedirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool recurseInNewSet=true)
Recursively replace all servers with the new servers in newSet.
const TNamed * namePtr() const
De-duplicated pointer to this object's name.
Definition RooAbsArg.h:561
void setShapeDirty()
Notify that a shape-like property (e.g. binning) has changed.
Definition RooAbsArg.h:493
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
bool isValueDirtyAndClear() const
Definition RooAbsArg.h:434
bool _fast
Definition RooAbsArg.h:715
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
RooFit::OwningPtr< RooArgSet > getComponents() const
Create a RooArgSet with all components (branch nodes) of the expression tree headed by this object.
bool dependsOnValue(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=nullptr) const
Check whether this object depends on values from an element in the serverList.
Definition RooAbsArg.h:106
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition RooAbsArg.h:488
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Implement multi-line detailed printing.
virtual RooAbsArg * cloneTree(const char *newname=nullptr) const
Clone tree expression of objects.
TString cleanBranchName() const
Construct a mangled name from the actual name that is free of any math symbols that might be interpre...
Int_t numProxies() const
Return the number of registered proxies.
static bool _inhibitDirty
Definition RooAbsArg.h:694
void setAttribute(const Text_t *name, bool value=true)
Set (default) or clear a named boolean attribute of this object.
void setProxyNormSet(const RooArgSet *nset)
Forward a change in the cached normalization argset to all the registered proxies.
void branchNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=nullptr, bool recurseNonDerived=false) const
Fill supplied list with all branch nodes of the arg tree starting with ourself as top node.
RooAbsProxy * getProxy(Int_t index) const
Return the nth proxy from the proxy list.
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition RooAbsArg.h:89
RefCountList_t _serverList
Definition RooAbsArg.h:632
void leafNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=nullptr, bool recurseNonDerived=false) const
Fill supplied list with all leaf nodes of the arg tree, starting with ourself as top node.
virtual bool isFundamental() const
Is this object a fundamental type that can be added to a dataset? Fundamental-type subclasses overrid...
Definition RooAbsArg.h:249
virtual bool redirectServersHook(const RooAbsCollection &newServerList, bool mustReplaceAll, bool nameChange, bool isRecursiveStep)
Function that is called at the end of redirectServers().
virtual bool checkObservables(const RooArgSet *nset) const
Overloadable function in which derived classes can implement consistency checks of the variables.
void treeNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=nullptr, bool doBranch=true, bool doLeaf=true, bool valueOnly=false, bool recurseNonDerived=false) const
Fill supplied list with nodes of the arg tree, following all server links, starting with ourself as t...
Abstract base class for RooRealVar binning definitions.
virtual bool isParameterized() const
Interface function.
virtual RooAbsReal * highBoundFunc() const
Return pointer to RooAbsReal parameterized upper bound, if any.
virtual RooAbsReal * lowBoundFunc() const
Return pointer to RooAbsReal parameterized lower bound, if any.
Abstract base class for objects that represent a discrete value that can be set from the outside,...
virtual bool setIndex(value_type index, bool printError=true)=0
Change category state by specifying the index code of the desired state.
A space to attach TBranches.
virtual value_type getCurrentIndex() const
Return index number of current state.
bool isSignType(bool mustHaveZero=false) const
Determine if category has 2 or 3 states with index values -1,0,1.
Abstract container object that can hold multiple RooAbsArg objects.
RooFit::UniqueId< RooAbsCollection > const & uniqueId() const
Returns a unique ID that is different for every instantiated RooAbsCollection.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
const char * GetName() const override
Returns name of object.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
bool overlaps(Iterator_t otherCollBegin, Iterator_t otherCollEnd) const
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
Storage_t::size_type size() const
RooAbsArg * first() const
virtual bool replace(const RooAbsArg &var1, const RooAbsArg &var2)
Replace var1 with var2 and return true for success.
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
RooAbsCollection * selectByName(const char *nameList, bool verbose=false) const
Create a subset of the current collection, consisting only of those elements with names matching the ...
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
void sort(bool reverse=false)
Sort collection using std::sort and name comparison.
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
std::string contentsString() const
Return comma separated list of contained object names as STL string.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
virtual double getMaxLimit(UInt_t dimension) const =0
virtual double operator()(const double xvector[]) const =0
virtual double getMinLimit(UInt_t dimension) const =0
Abstract base class for objects that are lvalues, i.e.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
static TClass * Class()
@ CanNotBeExtended
Definition RooAbsPdf.h:213
const char * normRange() const
Definition RooAbsPdf.h:251
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:57
virtual ExtendMode extendMode() const
Returns ability of PDF to provide extended likelihood terms.
Definition RooAbsPdf.h:217
Abstract interface for proxy classes.
Definition RooAbsProxy.h:37
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
TH1 * createHistogram(const char *name, 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
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
void setServerValues(const char *tmp)
Definition RooAbsReal.h:316
void setMessage(const char *tmp)
Definition RooAbsReal.h:315
Ref(RooAbsReal &ref)
Definition RooAbsReal.h:70
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, double scaleFactor, bool