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