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 RooAbsReal is the common 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 "RooChi2Var.h"
47#include "RooCmdConfig.h"
48#include "RooConstVar.h"
49#include "RooCurve.h"
50#include "RooCustomizer.h"
51#include "RooDataHist.h"
52#include "RooDataSet.h"
53#include "RooDerivative.h"
54#include "RooFirstMoment.h"
56#include "RooFit/Evaluator.h"
57#include "RooFitResult.h"
58#include "RooFormulaVar.h"
59#include "RooFunctor.h"
60#include "RooGlobalFunc.h"
61#include "RooHelpers.h"
62#include "RooHist.h"
63#include "RooMoment.h"
64#include "RooMsgService.h"
65#include "RooNumIntConfig.h"
66#include "RooNumRunningInt.h"
67#include "RooParamBinning.h"
68#include "RooPlot.h"
69#include "RooProduct.h"
70#include "RooProfileLL.h"
71#include "RooRealBinding.h"
72#include "RooRealIntegral.h"
73#include "RooRealVar.h"
74#include "RooSecondMoment.h"
75#include "RooVectorDataStore.h"
76#include "RooXYChi2Var.h"
77#include "TreeReadBuffer.h"
78#include "ValueChecking.h"
79
80#include "ROOT/StringUtils.hxx"
81#include "Compression.h"
82#include "Math/IFunction.h"
83#include "TMath.h"
84#include "TObjString.h"
85#include "TTree.h"
86#include "TH1.h"
87#include "TH2.h"
88#include "TH3.h"
89#include "TBranch.h"
90#include "TLeaf.h"
91#include "TAttLine.h"
92#include "TF1.h"
93#include "TF2.h"
94#include "TF3.h"
95#include "TMatrixD.h"
96#include "TVector.h"
97#include "strlcpy.h"
98#ifndef NDEBUG
99#include <TSystem.h> // To print stack traces when caching errors are detected
100#endif
101
102#include <iomanip>
103#include <iostream>
104#include <limits>
105#include <sstream>
106#include <sys/types.h>
107
108namespace {
109
110// Internal helper RooAbsFunc that evaluates the scaled data-weighted average of
111// given RooAbsReal as a function of a single variable using the RooFit::Evaluator.
112class ScaledDataWeightedAverage : public RooAbsFunc {
113public:
114 ScaledDataWeightedAverage(RooAbsReal const &arg, RooAbsData const &data, double scaleFactor, RooAbsRealLValue &var)
115 : RooAbsFunc{1}, _var{var}, _dataWeights{data.getWeightBatch(0, data.numEntries())}, _scaleFactor{scaleFactor}
116 {
117 _arg = RooFit::Detail::compileForNormSet(arg, *data.get());
118 _arg->recursiveRedirectServers(RooArgList{var});
119 _evaluator = std::make_unique<RooFit::Evaluator>(*_arg);
120 std::stack<std::vector<double>>{}.swap(_vectorBuffers);
121 auto dataSpans = RooFit::BatchModeDataHelpers::getDataSpans(data, "", nullptr, /*skipZeroWeights=*/false,
122 /*takeGlobalObservablesFromData=*/true,
123 _vectorBuffers);
124 for (auto const& item : dataSpans) {
125 _evaluator->setInput(item.first->GetName(), item.second, false);
126 }
127 }
128
129 double operator()(const double xvector[]) const override
130 {
131 double oldVal = _var.getVal();
132 _var.setVal(xvector[0]);
133
134 double out = 0.0;
135 std::span<const double> pdfValues = _evaluator->run();
136 if (_dataWeights.empty()) {
137 out = std::accumulate(pdfValues.begin(), pdfValues.end(), 0.0) / pdfValues.size();
138 } else {
139 double weightsSum = 0.0;
140 for (std::size_t i = 0; i < pdfValues.size(); ++i) {
141 out += pdfValues[i] * _dataWeights[i];
142 weightsSum += _dataWeights[i];
143 }
144 out /= weightsSum;
145 }
146 out *= _scaleFactor;
147
148 _var.setVal(oldVal);
149 return out;
150 }
151 double getMinLimit(UInt_t /*dimension*/) const override { return _var.getMin(); }
152 double getMaxLimit(UInt_t /*dimension*/) const override { return _var.getMax(); }
153
154private:
155 RooAbsRealLValue &_var;
156 std::unique_ptr<RooAbsReal> _arg;
157 std::span<const double> _dataWeights;
158 double _scaleFactor;
159 std::unique_ptr<RooFit::Evaluator> _evaluator;
160 std::stack<std::vector<double>> _vectorBuffers;
161};
162
163} // namespace
164
165using namespace std ;
166
168
170bool RooAbsReal::_hideOffset = true ;
171
172void RooAbsReal::setHideOffset(bool flag) { _hideOffset = flag ; }
174
177std::map<const RooAbsArg*,std::pair<std::string,std::list<RooAbsReal::EvalError> > > RooAbsReal::_evalErrorList ;
178
179
180////////////////////////////////////////////////////////////////////////////////
181/// coverity[UNINIT_CTOR]
182/// Default constructor
183
185
186
187////////////////////////////////////////////////////////////////////////////////
188/// Constructor with unit label
189
190RooAbsReal::RooAbsReal(const char *name, const char *title, const char *unit) :
191 RooAbsArg(name,title), _unit(unit)
192{
193 setValueDirty() ;
194 setShapeDirty() ;
195}
196
197
198////////////////////////////////////////////////////////////////////////////////
199/// Constructor with plot range and unit label
200
201RooAbsReal::RooAbsReal(const char *name, const char *title, double inMinVal,
202 double inMaxVal, const char *unit) :
203 RooAbsArg(name,title), _plotMin(inMinVal), _plotMax(inMaxVal), _unit(unit)
204{
205 setValueDirty() ;
206 setShapeDirty() ;
207}
208
209
210////////////////////////////////////////////////////////////////////////////////
211/// Copy constructor
212RooAbsReal::RooAbsReal(const RooAbsReal& other, const char* name) :
213 RooAbsArg(other,name), _plotMin(other._plotMin), _plotMax(other._plotMax),
214 _plotBins(other._plotBins), _value(other._value), _unit(other._unit), _label(other._label),
215 _forceNumInt(other._forceNumInt), _selectComp(other._selectComp)
216{
217 if (other._specIntegratorConfig) {
218 _specIntegratorConfig = std::make_unique<RooNumIntConfig>(*other._specIntegratorConfig) ;
219 }
220}
221
222
223////////////////////////////////////////////////////////////////////////////////
224/// Destructor
225
227
228
229////////////////////////////////////////////////////////////////////////////////
230/// Equality operator comparing to a double
231
233{
234 return (getVal()==value) ;
235}
236
237
238
239////////////////////////////////////////////////////////////////////////////////
240/// Equality operator when comparing to another RooAbsArg.
241/// Only functional when the other arg is a RooAbsReal
242
243bool RooAbsReal::operator==(const RooAbsArg& other) const
244{
245 const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
246 return otherReal ? operator==(otherReal->getVal()) : false ;
247}
248
249
250////////////////////////////////////////////////////////////////////////////////
251
252bool RooAbsReal::isIdentical(const RooAbsArg& other, bool assumeSameType) const
253{
254 if (!assumeSameType) {
255 const RooAbsReal* otherReal = dynamic_cast<const RooAbsReal*>(&other) ;
256 return otherReal ? operator==(otherReal->getVal()) : false ;
257 } else {
258 return getVal() == static_cast<const RooAbsReal&>(other).getVal();
259 }
260}
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// Return this variable's title string. If appendUnit is true and
265/// this variable has units, also append a string " (<unit>)".
266
267TString RooAbsReal::getTitle(bool appendUnit) const
268{
269 TString title(GetTitle());
270 if(appendUnit && 0 != strlen(getUnit())) {
271 title.Append(" (");
272 title.Append(getUnit());
273 title.Append(")");
274 }
275 return title;
276}
277
278
279
280////////////////////////////////////////////////////////////////////////////////
281/// Return value of object. If the cache is clean, return the
282/// cached value, otherwise recalculate on the fly and refill
283/// the cache
284
285double RooAbsReal::getValV(const RooArgSet* nset) const
286{
287 if (nset && nset->uniqueId().value() != _lastNormSetId) {
288 const_cast<RooAbsReal*>(this)->setProxyNormSet(nset);
289 _lastNormSetId = nset->uniqueId().value();
290 }
291
292 if (isValueDirtyAndClear()) {
293 _value = traceEval(nullptr) ;
294 // clearValueDirty() ;
295 }
296 // cout << "RooAbsReal::getValV(" << GetName() << ") writing _value = " << _value << std::endl ;
297
298 return hideOffset() ? _value + offset() : _value;
299}
300
301
302////////////////////////////////////////////////////////////////////////////////
303
305{
306 return _evalErrorList.size() ;
307}
308
309
310////////////////////////////////////////////////////////////////////////////////
311
313{
314 return _evalErrorList.begin() ;
315}
316
317
318////////////////////////////////////////////////////////////////////////////////
319/// Calculate current value of object, with error tracing wrapper
320
321double RooAbsReal::traceEval(const RooArgSet* /*nset*/) const
322{
323 double value = evaluate() ;
324
325 if (TMath::IsNaN(value)) {
326 logEvalError("function value is NAN") ;
327 }
328
329 //cxcoutD(Tracing) << "RooAbsReal::getValF(" << GetName() << ") operMode = " << _operMode << " recalculated, new value = " << value << std::endl ;
330
331 //Standard tracing code goes here
332 if (!isValidReal(value)) {
333 coutW(Tracing) << "RooAbsReal::traceEval(" << GetName()
334 << "): validation failed: " << value << std::endl ;
335 }
336
337 //Call optional subclass tracing code
338 // traceEvalHook(value) ;
339
340 return value ;
341}
342
343
344
345////////////////////////////////////////////////////////////////////////////////
346/// Variant of getAnalyticalIntegral that is also passed the normalization set
347/// that should be applied to the integrand of which the integral is requested.
348/// For certain operator p.d.f it is useful to overload this function rather
349/// than analyticalIntegralWN() as the additional normalization information
350/// may be useful in determining a more efficient decomposition of the
351/// requested integral.
352
354 const RooArgSet* /*normSet*/, const char* rangeName) const
355{
356 return _forceNumInt ? 0 : getAnalyticalIntegral(allDeps,analDeps,rangeName) ;
357}
358
359
360
361////////////////////////////////////////////////////////////////////////////////
362/// Interface function getAnalyticalIntergral advertises the
363/// analytical integrals that are supported. 'integSet'
364/// is the set of dependents for which integration is requested. The
365/// function should copy the subset of dependents it can analytically
366/// integrate to anaIntSet and return a unique identification code for
367/// this integration configuration. If no integration can be
368/// performed, zero should be returned.
369
370Int_t RooAbsReal::getAnalyticalIntegral(RooArgSet& /*integSet*/, RooArgSet& /*anaIntSet*/, const char* /*rangeName*/) const
371{
372 return 0 ;
373}
374
375
376
377////////////////////////////////////////////////////////////////////////////////
378/// Implements the actual analytical integral(s) advertised by
379/// getAnalyticalIntegral. This functions will only be called with
380/// codes returned by getAnalyticalIntegral, except code zero.
381
382double RooAbsReal::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
383{
384// cout << "RooAbsReal::analyticalIntegralWN(" << GetName() << ") code = " << code << " normSet = " << (normSet?*normSet:RooArgSet()) << std::endl ;
385 if (code==0) return getVal(normSet) ;
386 return analyticalIntegral(code,rangeName) ;
387}
388
389
390
391////////////////////////////////////////////////////////////////////////////////
392/// Implements the actual analytical integral(s) advertised by
393/// getAnalyticalIntegral. This functions will only be called with
394/// codes returned by getAnalyticalIntegral, except code zero.
395
396double RooAbsReal::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
397{
398 // By default no analytical integrals are implemented
399 coutF(Eval) << "RooAbsReal::analyticalIntegral(" << GetName() << ") code " << code << " not implemented" << std::endl ;
400 return 0 ;
401}
402
403
404
405////////////////////////////////////////////////////////////////////////////////
406/// Get the label associated with the variable
407
408const char *RooAbsReal::getPlotLabel() const
409{
410 return _label.IsNull() ? fName.Data() : _label.Data();
411}
412
413
414
415////////////////////////////////////////////////////////////////////////////////
416/// Set the label associated with this variable
417
418void RooAbsReal::setPlotLabel(const char *label)
419{
420 _label= label;
421}
422
423
424
425////////////////////////////////////////////////////////////////////////////////
426///Read object contents from stream (dummy for now)
427
428bool RooAbsReal::readFromStream(std::istream& /*is*/, bool /*compact*/, bool /*verbose*/)
429{
430 return false ;
431}
432
433
434
435////////////////////////////////////////////////////////////////////////////////
436///Write object contents to stream (dummy for now)
437
438void RooAbsReal::writeToStream(std::ostream& /*os*/, bool /*compact*/) const
439{
440}
441
442
443
444////////////////////////////////////////////////////////////////////////////////
445/// Print object value
446
447void RooAbsReal::printValue(std::ostream& os) const
448{
449 os << getVal() ;
450}
451
452
453
454////////////////////////////////////////////////////////////////////////////////
455/// Structure printing
456
457void RooAbsReal::printMultiline(std::ostream& os, Int_t contents, bool verbose, TString indent) const
458{
459 RooAbsArg::printMultiline(os,contents,verbose,indent) ;
460 os << indent << "--- RooAbsReal ---" << std::endl;
461 TString unit(_unit);
462 if(!unit.IsNull()) unit.Prepend(' ');
463 //os << indent << " Value = " << getVal() << unit << std::endl;
464 os << std::endl << indent << " Plot label is \"" << getPlotLabel() << "\"" << "\n";
465}
466
467
468////////////////////////////////////////////////////////////////////////////////
469/// Create a RooProfileLL object that eliminates all nuisance parameters in the
470/// present function. The nuisance parameters are defined as all parameters
471/// of the function except the stated paramsOfInterest
472
474{
475 // Construct name of profile object
476 auto name = std::string(GetName()) + "_Profile[";
477 bool first = true;
478 for (auto const& arg : paramsOfInterest) {
479 if (first) {
480 first = false ;
481 } else {
482 name.append(",") ;
483 }
484 name.append(arg->GetName()) ;
485 }
486 name.append("]") ;
487
488 // Create and return profile object
489 auto out = std::make_unique<RooProfileLL>(name.c_str(),(std::string("Profile of ") + GetTitle()).c_str(),*this,paramsOfInterest);
490 return RooFit::Detail::owningPtr(std::move(out));
491}
492
493
494
495
496
497
498////////////////////////////////////////////////////////////////////////////////
499/// Create an object that represents the integral of the function over one or more observables listed in `iset`.
500/// The actual integration calculation is only performed when the returned object is evaluated. The name
501/// of the integral object is automatically constructed from the name of the input function, the variables
502/// it integrates and the range integrates over.
503///
504/// \note The integral over a PDF is usually not normalised (*i.e.*, it is usually not
505/// 1 when integrating the PDF over the full range). In fact, this integral is used *to compute*
506/// the normalisation of each PDF. See the rf110 tutorial at https://root.cern.ch/doc/master/group__tutorial__roofit.html
507/// for details on PDF normalisation.
508///
509/// The following named arguments are accepted
510/// | | Effect on integral creation
511/// |--|-------------------------------
512/// | `NormSet(const RooArgSet&)` | Specify normalization set, mostly useful when working with PDFs
513/// | `NumIntConfig(const RooNumIntConfig&)` | Use given configuration for any numeric integration, if necessary
514/// | `Range(const char* name)` | Integrate only over given range. Multiple ranges may be specified by passing multiple Range() arguments
515
517 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
518 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
519{
520
521
522 // Define configuration for this method
523 RooCmdConfig pc("RooAbsReal::createIntegral(" + std::string(GetName()) + ")");
524 pc.defineString("rangeName","RangeWithName",0,"",true) ;
525 pc.defineSet("normSet","NormSet",0,nullptr) ;
526 pc.defineObject("numIntConfig","NumIntConfig",0,nullptr) ;
527
528 // Process & check varargs
529 pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
530 if (!pc.ok(true)) {
531 return nullptr;
532 }
533
534 // Extract values from named arguments
535 const char* rangeName = pc.getString("rangeName",nullptr,true) ;
536 const RooArgSet* nset = pc.getSet("normSet",nullptr);
537 const RooNumIntConfig* cfg = static_cast<const RooNumIntConfig*>(pc.getObject("numIntConfig",nullptr)) ;
538
539 return createIntegral(iset,nset,cfg,rangeName) ;
540}
541
542
543
544
545
546////////////////////////////////////////////////////////////////////////////////
547/// Create an object that represents the integral of the function over one or more observables listed in iset.
548/// The actual integration calculation is only performed when the return object is evaluated. The name
549/// of the integral object is automatically constructed from the name of the input function, the variables
550/// it integrates and the range integrates over. If nset is specified the integrand is request
551/// to be normalized over nset (only meaningful when the integrand is a pdf). If rangename is specified
552/// the integral is performed over the named range, otherwise it is performed over the domain of each
553/// integrated observable. If cfg is specified it will be used to configure any numeric integration
554/// aspect of the integral. It will not force the integral to be performed numerically, which is
555/// decided automatically by RooRealIntegral.
556
558 const RooNumIntConfig* cfg, const char* rangeName) const
559{
560 if (!rangeName || strchr(rangeName,',')==nullptr) {
561 // Simple case: integral over full range or single limited range
562 return createIntObj(iset,nset,cfg,rangeName);
563 }
564
565 // Integral over multiple ranges
566 std::vector<std::string> tokens = ROOT::Split(rangeName, ",");
567
568 if(RooHelpers::checkIfRangesOverlap(iset, tokens)) {
569 std::stringstream errMsg;
570 errMsg << GetName() << " : integrating with respect to the variables " << iset << " on the ranges \"" << rangeName
571 << "\" is not possible because the ranges are overlapping";
572 const std::string errMsgString = errMsg.str();
573 coutE(Integration) << errMsgString << std::endl;
574 throw std::invalid_argument(errMsgString);
575 }
576
577 RooArgSet components ;
578 for (const std::string& token : tokens) {
579 components.addOwned(std::unique_ptr<RooAbsReal>{createIntObj(iset,nset,cfg, token.c_str())});
580 }
581
582 const std::string title = std::string("Integral of ") + GetTitle();
583 const std::string fullName = std::string(GetName()) + integralNameSuffix(iset,nset,rangeName).Data();
584
585 auto out = std::make_unique<RooAddition>(fullName.c_str(), title.c_str(), components);
586 out->addOwnedComponents(std::move(components));
587 return RooFit::Detail::owningPtr<RooAbsReal>(std::move(out));
588}
589
590
591
592////////////////////////////////////////////////////////////////////////////////
593/// Internal utility function for createIntegral() that creates the actual integral object.
595 const RooNumIntConfig* cfg, const char* rangeName) const
596{
597 // Make internal use copies of iset and nset
598 RooArgSet iset(iset2) ;
599 const RooArgSet* nset = nset2 ;
600
601
602 // Initialize local variables perparing for recursive loop
603 bool error = false ;
604 const RooAbsReal* integrand = this ;
605 std::unique_ptr<RooAbsReal> integral;
606
607 // Handle trivial case of no integration here explicitly
608 if (iset.empty()) {
609
610 const std::string title = std::string("Integral of ") + GetTitle();
611 const std::string name = std::string(GetName()) + integralNameSuffix(iset,nset,rangeName).Data();
612
613 auto out = std::make_unique<RooRealIntegral>(name.c_str(), title.c_str(), *this, iset, nset, cfg, rangeName);
614 return RooFit::Detail::owningPtr<RooAbsReal>(std::move(out));
615 }
616
617 // Process integration over remaining integration variables
618 while(!iset.empty()) {
619
620
621 // Find largest set of observables that can be integrated in one go
622 RooArgSet innerSet ;
623 findInnerMostIntegration(iset,innerSet,rangeName) ;
624
625 // If largest set of observables that can be integrated is empty set, problem was ill defined
626 // Postpone error messaging and handling to end of function, exit loop here
627 if (innerSet.empty()) {
628 error = true ;
629 break ;
630 }
631
632 // Prepare name and title of integral to be created
633 const std::string title = std::string("Integral of ") + integrand->GetTitle();
634 const std::string name = std::string(integrand->GetName()) + integrand->integralNameSuffix(innerSet,nset,rangeName).Data();
635
636 std::unique_ptr<RooAbsReal> innerIntegral = std::move(integral);
637
638 // Construct innermost integral
639 integral = std::make_unique<RooRealIntegral>(name.c_str(),title.c_str(),*integrand,innerSet,nset,cfg,rangeName);
640
641 // Integral of integral takes ownership of innermost integral
642 if (innerIntegral) {
643 integral->addOwnedComponents(std::move(innerIntegral));
644 }
645
646 // Remove already integrated observables from to-do list
647 iset.remove(innerSet) ;
648
649 // Send info message on recursion if needed
650 if (integrand == this && !iset.empty()) {
651 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 ;
652 }
653
654 // Prepare for recursion, next integral should integrate last integrand
655 integrand = integral.get();
656
657
658 // Only need normalization set in innermost integration
659 nset = nullptr;
660 }
661
662 if (error) {
663 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;
664 return nullptr;
665 }
666
667
668 // After-burner: apply interpolating cache on (numeric) integral if requested by user
669 const char* cacheParamsStr = getStringAttribute("CACHEPARAMINT") ;
670 if (cacheParamsStr && strlen(cacheParamsStr)) {
671
672 std::unique_ptr<RooArgSet> intParams{integral->getVariables()};
673
674 RooArgSet cacheParams = RooHelpers::selectFromArgSet(*intParams, cacheParamsStr);
675
676 if (!cacheParams.empty()) {
677 cxcoutD(Caching) << "RooAbsReal::createIntObj(" << GetName() << ") INFO: constructing " << cacheParams.size()
678 << "-dim value cache for integral over " << iset2 << " as a function of " << cacheParams << " in range " << (rangeName?rangeName:"<none>") << std::endl ;
679 std::string name = Form("%s_CACHE_[%s]",integral->GetName(),cacheParams.contentsString().c_str()) ;
680 auto cachedIntegral = std::make_unique<RooCachedReal>(name.c_str(),name.c_str(),*integral,cacheParams);
681 cachedIntegral->setInterpolationOrder(2) ;
682 cachedIntegral->addOwnedComponents(std::move(integral));
683 cachedIntegral->setCacheSource(true) ;
684 if (integral->operMode()==ADirty) {
685 cachedIntegral->setOperMode(ADirty) ;
686 }
687 //cachedIntegral->disableCache(true) ;
688 return RooFit::Detail::owningPtr<RooAbsReal>(std::move(cachedIntegral));
689 }
690 }
691
692 return RooFit::Detail::owningPtr(std::move(integral));
693}
694
695
696
697////////////////////////////////////////////////////////////////////////////////
698/// Utility function for createIntObj() that aids in the construct of recursive integrals
699/// over functions with multiple observables with parameterized ranges. This function
700/// finds in a given set allObs over which integration is requested the largeset subset
701/// of observables that can be integrated simultaneously. This subset consists of
702/// observables with fixed ranges and observables with parameterized ranges whose
703/// parameterization does not depend on any observable that is also integrated.
704
705void RooAbsReal::findInnerMostIntegration(const RooArgSet& allObs, RooArgSet& innerObs, const char* rangeName) const
706{
707 // Make lists of
708 // a) integrated observables with fixed ranges,
709 // b) integrated observables with parameterized ranges depending on other integrated observables
710 // c) integrated observables used in definition of any parameterized ranges of integrated observables
711 RooArgSet obsWithFixedRange(allObs) ;
712 RooArgSet obsWithParamRange ;
713 RooArgSet obsServingAsRangeParams ;
714
715 // Loop over all integrated observables
716 for (const auto aarg : allObs) {
717 // Check if observable is real-valued lvalue
718 if (auto arglv = dynamic_cast<RooAbsRealLValue*>(aarg)) {
719
720 // Check if range is parameterized
721 RooAbsBinning& binning = arglv->getBinning(rangeName,false,true) ;
722 if (binning.isParameterized()) {
723 RooArgSet loBoundObs;
724 RooArgSet hiBoundObs;
725 binning.lowBoundFunc()->getObservables(&allObs, loBoundObs) ;
726 binning.highBoundFunc()->getObservables(&allObs, hiBoundObs) ;
727
728 // Check if range parameterization depends on other integrated observables
729 if (loBoundObs.overlaps(allObs) || hiBoundObs.overlaps(allObs)) {
730 obsWithParamRange.add(*aarg) ;
731 obsWithFixedRange.remove(*aarg) ;
732 obsServingAsRangeParams.add(loBoundObs,false) ;
733 obsServingAsRangeParams.add(hiBoundObs,false) ;
734 }
735 }
736 }
737 }
738
739 // Make list of fixed-range observables that are _not_ involved in the parameterization of ranges of other observables
740 RooArgSet obsWithFixedRangeNP(obsWithFixedRange) ;
741 obsWithFixedRangeNP.remove(obsServingAsRangeParams) ;
742
743 // Make list of param-range observables that are _not_ involved in the parameterization of ranges of other observables
744 RooArgSet obsWithParamRangeNP(obsWithParamRange) ;
745 obsWithParamRangeNP.remove(obsServingAsRangeParams) ;
746
747 // Construct inner-most integration: over observables (with fixed or param range) not used in any other param range definitions
748 innerObs.removeAll() ;
749 innerObs.add(obsWithFixedRangeNP) ;
750 innerObs.add(obsWithParamRangeNP) ;
751
752}
753
754
755////////////////////////////////////////////////////////////////////////////////
756/// Construct string with unique suffix name to give to integral object that encodes
757/// integrated observables, normalization observables and the integration range name
758
759TString RooAbsReal::integralNameSuffix(const RooArgSet& iset, const RooArgSet* nset, const char* rangeName, bool omitEmpty) const
760{
761 TString name ;
762 if (!iset.empty()) {
763
764 RooArgSet isetTmp(iset) ;
765 isetTmp.sort() ;
766
767 name.Append("_Int[") ;
768 bool first(true) ;
769 for(RooAbsArg * arg : isetTmp) {
770 if (first) {
771 first=false ;
772 } else {
773 name.Append(",") ;
774 }
775 name.Append(arg->GetName()) ;
776 }
777 if (rangeName) {
778 name.Append("|") ;
779 name.Append(rangeName) ;
780 }
781 name.Append("]");
782 } else if (!omitEmpty) {
783 name.Append("_Int[]") ;
784 }
785
786 if (nset && !nset->empty()) {
787
788 RooArgSet nsetTmp(*nset) ;
789 nsetTmp.sort() ;
790
791 name.Append("_Norm[") ;
792 bool first(true);
793 for(RooAbsArg * arg : nsetTmp) {
794 if (first) {
795 first=false ;
796 } else {
797 name.Append(",") ;
798 }
799 name.Append(arg->GetName()) ;
800 }
801 const RooAbsPdf* thisPdf = dynamic_cast<const RooAbsPdf*>(this) ;
802 if (thisPdf && thisPdf->normRange()) {
803 name.Append("|") ;
804 name.Append(thisPdf->normRange()) ;
805 }
806 name.Append("]") ;
807 }
808
809 return name ;
810}
811
812
813
814////////////////////////////////////////////////////////////////////////////////
815/// Utility function for plotOn() that creates a projection of a function or p.d.f
816/// to be plotted on a RooPlot.
817/// \ref createPlotProjAnchor "createPlotProjection()"
818
819const RooAbsReal* RooAbsReal::createPlotProjection(const RooArgSet& depVars, const RooArgSet& projVars,
820 RooArgSet*& cloneSet) const
821{
822 return createPlotProjection(depVars,&projVars,cloneSet) ;
823}
824
825
826////////////////////////////////////////////////////////////////////////////////
827/// Utility function for plotOn() that creates a projection of a function or p.d.f
828/// to be plotted on a RooPlot.
829/// \anchor createPlotProjAnchor
830///
831/// Create a new object \f$ G \f$ that represents the normalized projection:
832/// \f[
833/// G[x,p] = \frac{\int F[x,y,p] \; \mathrm{d}\{y\}}
834/// {\int F[x,y,p] \; \mathrm{d}\{x\} \, \mathrm{d}\{y\}}
835/// \f]
836/// where \f$ F[x,y,p] \f$ is the function we represent, and
837/// \f$ \{ p \} \f$ are the remaining variables ("parameters").
838///
839/// \param[in] dependentVars Dependent variables over which to normalise, \f$ \{x\} \f$.
840/// \param[in] projectedVars Variables to project out, \f$ \{ y \} \f$.
841/// \param[out] cloneSet Will be set to a RooArgSet*, which will contain a clone of *this plus its projection integral object.
842/// The latter will also be returned. The caller takes ownership of this set.
843/// \param[in] rangeName Optional range for projection integrals
844/// \param[in] condObs Conditional observables, which are not integrated for normalisation, even if they
845/// are in `dependentVars` or `projectedVars`.
846/// \return A pointer to the newly created object, or zero in case of an
847/// error. The caller is responsible for deleting the `cloneSet` (which includes the returned projection object).
848const RooAbsReal *RooAbsReal::createPlotProjection(const RooArgSet &dependentVars, const RooArgSet *projectedVars,
849 RooArgSet *&cloneSet, const char* rangeName, const RooArgSet* condObs) const
850{
851 // Get the set of our leaf nodes
852 RooArgSet leafNodes;
853 RooArgSet treeNodes;
854 leafNodeServerList(&leafNodes,this);
855 treeNodeServerList(&treeNodes,this) ;
856
857
858 // Check that the dependents are all fundamental. Filter out any that we
859 // do not depend on, and make substitutions by name in our leaf list.
860 // Check for overlaps with the projection variables.
861 for (const auto arg : dependentVars) {
862 if(!arg->isFundamental() && !dynamic_cast<const RooAbsLValue*>(arg)) {
863 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: variable \"" << arg->GetName()
864 << "\" of wrong type: " << arg->ClassName() << std::endl;
865 return nullptr;
866 }
867
868 RooAbsArg *found= treeNodes.find(arg->GetName());
869 if(!found) {
870 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
871 << "\" is not a dependent and will be ignored." << std::endl;
872 continue;
873 }
874 if(found != arg) {
875 if (leafNodes.find(found->GetName())) {
876 leafNodes.replace(*found,*arg);
877 } else {
878 leafNodes.add(*arg) ;
879
880 // Remove any dependents of found, replace by dependents of LV node
881 RooArgSet lvDep;
882 arg->getObservables(&leafNodes, lvDep);
883 for (const auto lvs : lvDep) {
884 RooAbsArg* tmp = leafNodes.find(lvs->GetName()) ;
885 if (tmp) {
886 leafNodes.remove(*tmp) ;
887 leafNodes.add(*lvs) ;
888 }
889 }
890 }
891 }
892
893 // check if this arg is also in the projection set
894 if(nullptr != projectedVars && projectedVars->find(arg->GetName())) {
895 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: \"" << arg->GetName()
896 << "\" cannot be both a dependent and a projected variable." << std::endl;
897 return nullptr;
898 }
899 }
900
901 // Remove the projected variables from the list of leaf nodes, if necessary.
902 if(nullptr != projectedVars) leafNodes.remove(*projectedVars,true);
903
904 // Make a deep-clone of ourself so later operations do not disturb our original state
905 cloneSet = new RooArgSet;
906 if (RooArgSet(*this).snapshot(*cloneSet, true)) {
907 coutE(Plotting) << "RooAbsPdf::createPlotProjection(" << GetName() << ") Couldn't deep-clone PDF, abort," << std::endl ;
908 return nullptr ;
909 }
910 RooAbsReal *theClone= (RooAbsReal*)cloneSet->find(GetName());
911
912 // The remaining entries in our list of leaf nodes are the external
913 // dependents (x) and parameters (p) of the projection. Patch them back
914 // into the theClone. This orphans the nodes they replace, but the orphans
915 // are still in the cloneList and so will be cleaned up eventually.
916 //cout << "redirection leafNodes : " ; leafNodes.Print("1") ;
917
918 std::unique_ptr<RooArgSet> plotLeafNodes{static_cast<RooArgSet*>(leafNodes.selectCommon(dependentVars))};
919 theClone->recursiveRedirectServers(*plotLeafNodes,false,false,false);
920
921 // Create the set of normalization variables to use in the projection integrand
922 RooArgSet normSet(dependentVars);
923 if(nullptr != projectedVars) normSet.add(*projectedVars);
924 if(nullptr != condObs) {
925 normSet.remove(*condObs,true,true) ;
926 }
927
928 // Try to create a valid projection integral. If no variables are to be projected,
929 // create a null projection anyway to bind our normalization over the dependents
930 // consistently with the way they would be bound with a non-trivial projection.
931 RooArgSet empty;
932 if(nullptr == projectedVars) projectedVars= &empty;
933
934 TString name = GetName() ;
935 name += integralNameSuffix(*projectedVars,&normSet,rangeName,true) ;
936
937 TString title(GetTitle());
938 title.Prepend("Projection of ");
939
940
941 std::unique_ptr<RooAbsReal> projected{theClone->createIntegral(*projectedVars,normSet,rangeName)};
942
943 if(nullptr == projected || !projected->isValid()) {
944 coutE(Plotting) << ClassName() << "::" << GetName() << ":createPlotProjection: cannot integrate out ";
945 projectedVars->printStream(std::cout,kName|kArgs,kSingleLine);
946 return nullptr;
947 }
948
949 if(projected->InheritsFrom(RooRealIntegral::Class())){
950 static_cast<RooRealIntegral&>(*projected).setAllowComponentSelection(true);
951 }
952
953 projected->SetName(name.Data()) ;
954 projected->SetTitle(title.Data()) ;
955
956 // Add the projection integral to the cloneSet so that it eventually gets cleaned up by the caller.
957 RooAbsReal *projectedPtr = projected.get();
958 cloneSet->addOwned(std::move(projected));
959
960 // return a const pointer to remind the caller that they do not delete the returned object
961 // directly (it is contained in the cloneSet instead).
962 return projectedPtr;
963}
964
965
966
967////////////////////////////////////////////////////////////////////////////////
968/// Fill the ROOT histogram 'hist' with values sampled from this
969/// function at the bin centers. Our value is calculated by first
970/// integrating out any variables in projectedVars and then scaling
971/// the result by scaleFactor. Returns a pointer to the input
972/// histogram, or zero in case of an error. The input histogram can
973/// be any TH1 subclass, and therefore of arbitrary
974/// dimension. Variables are matched with the (x,y,...) dimensions of
975/// the input histogram according to the order in which they appear
976/// in the input plotVars list. If scaleForDensity is true the
977/// histogram is filled with a the functions density rather than
978/// the functions value (i.e. the value at the bin center is multiplied
979/// with bin volume)
980
982 double scaleFactor, const RooArgSet *projectedVars, bool scaleForDensity,
983 const RooArgSet* condObs, bool setError) const
984{
985 // Do we have a valid histogram to use?
986 if(nullptr == hist) {
987 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: no valid histogram to fill" << std::endl;
988 return nullptr;
989 }
990
991 // Check that the number of plotVars matches the input histogram's dimension
992 Int_t hdim= hist->GetDimension();
993 if(hdim != plotVars.getSize()) {
994 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: plotVars has the wrong dimension" << std::endl;
995 return nullptr;
996 }
997
998
999 // Check that the plot variables are all actually RooRealVars and print a warning if we do not
1000 // explicitly depend on one of them. Fill a set (not list!) of cloned plot variables.
1001 RooArgSet plotClones;
1002 for(Int_t index= 0; index < plotVars.getSize(); index++) {
1003 const RooAbsArg *var= plotVars.at(index);
1004 const RooRealVar *realVar= dynamic_cast<const RooRealVar*>(var);
1005 if(nullptr == realVar) {
1006 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot plot variable \"" << var->GetName()
1007 << "\" of type " << var->ClassName() << std::endl;
1008 return nullptr;
1009 }
1010 if(!this->dependsOn(*realVar)) {
1011 coutE(InputArguments) << ClassName() << "::" << GetName()
1012 << ":fillHistogram: WARNING: variable is not an explicit dependent: " << realVar->GetName() << std::endl;
1013 }
1014 plotClones.addClone(*realVar,true); // do not complain about duplicates
1015 }
1016
1017 // Reconnect all plotClones to each other, imported when plotting N-dim integrals with entangled parameterized ranges
1018 for(RooAbsArg * pc : plotClones) {
1019 pc->recursiveRedirectServers(plotClones,false,false,true) ;
1020 }
1021
1022 // Call checkObservables
1023 RooArgSet allDeps(plotClones) ;
1024 if (projectedVars) {
1025 allDeps.add(*projectedVars) ;
1026 }
1027 if (checkObservables(&allDeps)) {
1028 coutE(InputArguments) << "RooAbsReal::fillHistogram(" << GetName() << ") error in checkObservables, abort" << std::endl ;
1029 return hist ;
1030 }
1031
1032 // Create a standalone projection object to use for calculating bin contents
1033 RooArgSet *cloneSet = nullptr;
1034 const RooAbsReal *projected= createPlotProjection(plotClones,projectedVars,cloneSet,nullptr,condObs);
1035
1036 cxcoutD(Plotting) << "RooAbsReal::fillHistogram(" << GetName() << ") plot projection object is " << projected->GetName() << std::endl ;
1037
1038 // Prepare to loop over the histogram bins
1039 Int_t xbins(0),ybins(1),zbins(1);
1040 RooRealVar *xvar = nullptr;
1041 RooRealVar *yvar = nullptr;
1042 RooRealVar *zvar = nullptr;
1043 TAxis *xaxis = nullptr;
1044 TAxis *yaxis = nullptr;
1045 TAxis *zaxis = nullptr;
1046 switch(hdim) {
1047 case 3:
1048 zbins= hist->GetNbinsZ();
1049 zvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(2)->GetName()));
1050 zaxis= hist->GetZaxis();
1051 assert(nullptr != zvar && nullptr != zaxis);
1052 if (scaleForDensity) {
1053 scaleFactor*= (zaxis->GetXmax() - zaxis->GetXmin())/zbins;
1054 }
1055 // fall through to next case...
1056 case 2:
1057 ybins= hist->GetNbinsY();
1058 yvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(1)->GetName()));
1059 yaxis= hist->GetYaxis();
1060 assert(nullptr != yvar && nullptr != yaxis);
1061 if (scaleForDensity) {
1062 scaleFactor*= (yaxis->GetXmax() - yaxis->GetXmin())/ybins;
1063 }
1064 // fall through to next case...
1065 case 1:
1066 xbins= hist->GetNbinsX();
1067 xvar= dynamic_cast<RooRealVar*>(plotClones.find(plotVars.at(0)->GetName()));
1068 xaxis= hist->GetXaxis();
1069 assert(nullptr != xvar && nullptr != xaxis);
1070 if (scaleForDensity) {
1071 scaleFactor*= (xaxis->GetXmax() - xaxis->GetXmin())/xbins;
1072 }
1073 break;
1074 default:
1075 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillHistogram: cannot fill histogram with "
1076 << hdim << " dimensions" << std::endl;
1077 break;
1078 }
1079
1080 // Loop over the input histogram's bins and fill each one with our projection's
1081 // value, calculated at the center.
1083 Int_t xbin(0),ybin(0),zbin(0);
1084 Int_t bins= xbins*ybins*zbins;
1085 for(Int_t bin= 0; bin < bins; bin++) {
1086 switch(hdim) {
1087 case 3:
1088 if(bin % (xbins*ybins) == 0) {
1089 zbin++;
1090 zvar->setVal(zaxis->GetBinCenter(zbin));
1091 }
1092 // fall through to next case...
1093 case 2:
1094 if(bin % xbins == 0) {
1095 ybin= (ybin%ybins) + 1;
1096 yvar->setVal(yaxis->GetBinCenter(ybin));
1097 }
1098 // fall through to next case...
1099 case 1:
1100 xbin= (xbin%xbins) + 1;
1101 xvar->setVal(xaxis->GetBinCenter(xbin));
1102 break;
1103 default:
1104 coutE(InputArguments) << "RooAbsReal::fillHistogram: Internal Error!" << std::endl;
1105 break;
1106 }
1107
1108 double result= scaleFactor*projected->getVal();
1109 if (RooAbsReal::numEvalErrors()>0) {
1110 coutW(Plotting) << "WARNING: Function evaluation error(s) at coordinates [x]=" << xvar->getVal() ;
1111 if (hdim==2) ccoutW(Plotting) << " [y]=" << yvar->getVal() ;
1112 if (hdim==3) ccoutW(Plotting) << " [z]=" << zvar->getVal() ;
1113 ccoutW(Plotting) << std::endl ;
1114 // RooAbsReal::printEvalErrors(ccoutW(Plotting),10) ;
1115 result = 0 ;
1116 }
1118
1119 hist->SetBinContent(hist->GetBin(xbin,ybin,zbin),result);
1120 if (setError) {
1121 hist->SetBinError(hist->GetBin(xbin,ybin,zbin),sqrt(result)) ;
1122 }
1123
1124 //cout << "bin " << bin << " -> (" << xbin << "," << ybin << "," << zbin << ") = " << result << std::endl;
1125 }
1127
1128 // cleanup
1129 delete cloneSet;
1130
1131 return hist;
1132}
1133
1134
1135
1136////////////////////////////////////////////////////////////////////////////////
1137/// Fill a RooDataHist with values sampled from this function at the
1138/// bin centers. If extendedMode is true, the p.d.f. values is multiplied
1139/// by the number of expected events in each bin
1140///
1141/// An optional scaling by a given scaleFactor can be performed.
1142/// Returns a pointer to the input RooDataHist, or zero
1143/// in case of an error.
1144///
1145/// If correctForBinSize is true the RooDataHist
1146/// is filled with the functions density (function value times the
1147/// bin volume) rather than function value.
1148///
1149/// If showProgress is true
1150/// a process indicator is printed on stdout in steps of one percent,
1151/// which is mostly useful for the sampling of expensive functions
1152/// such as likelihoods
1153
1154RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSet, double scaleFactor,
1155 bool correctForBinSize, bool showProgress) const
1156{
1157 // Do we have a valid histogram to use?
1158 if(nullptr == hist) {
1159 coutE(InputArguments) << ClassName() << "::" << GetName() << ":fillDataHist: no valid RooDataHist to fill" << std::endl;
1160 return nullptr;
1161 }
1162
1163 // Call checkObservables
1164 RooArgSet allDeps(*hist->get()) ;
1165 if (checkObservables(&allDeps)) {
1166 coutE(InputArguments) << "RooAbsReal::fillDataHist(" << GetName() << ") error in checkObservables, abort" << std::endl ;
1167 return hist ;
1168 }
1169
1170 // Make deep clone of self and attach to dataset observables
1171 //RooArgSet* origObs = getObservables(hist) ;
1172 RooArgSet cloneSet;
1173 RooArgSet(*this).snapshot(cloneSet, true);
1174 RooAbsReal* theClone = static_cast<RooAbsReal*>(cloneSet.find(GetName()));
1175 theClone->recursiveRedirectServers(*hist->get()) ;
1176 //const_cast<RooAbsReal*>(this)->recursiveRedirectServers(*hist->get()) ;
1177
1178 // Iterator over all bins of RooDataHist and fill weights
1179 Int_t onePct = hist->numEntries()/100 ;
1180 if (onePct==0) {
1181 onePct++ ;
1182 }
1183 for (Int_t i=0 ; i<hist->numEntries() ; i++) {
1184 if (showProgress && (i%onePct==0)) {
1185 ccoutP(Eval) << "." << std::flush ;
1186 }
1187 const RooArgSet* obs = hist->get(i) ;
1188 double binVal = theClone->getVal(normSet?normSet:obs)*scaleFactor ;
1189 if (correctForBinSize) {
1190 binVal*= hist->binVolume() ;
1191 }
1192 hist->set(i, binVal, 0.);
1193 }
1194
1195 return hist;
1196}
1197
1198
1199
1200
1201////////////////////////////////////////////////////////////////////////////////
1202/// Create and fill a ROOT histogram TH1, TH2 or TH3 with the values of this function for the variables with given names.
1203/// \param[in] varNameList List of variables to use for x, y, z axis, separated by ':'
1204/// \param[in] xbins Number of bins for first variable
1205/// \param[in] ybins Number of bins for second variable
1206/// \param[in] zbins Number of bins for third variable
1207/// \return TH1*, which is one of TH[1-3]. The histogram is owned by the caller.
1208///
1209/// For a greater degree of control use
1210/// RooAbsReal::createHistogram(const char *, const RooAbsRealLValue&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&, const RooCmdArg&) const
1211///
1212
1213TH1* RooAbsReal::createHistogram(RooStringView varNameList, Int_t xbins, Int_t ybins, Int_t zbins) const
1214{
1215 std::unique_ptr<RooArgSet> vars{getVariables()};
1216
1217 auto varNames = ROOT::Split(varNameList, ",:");
1218 std::vector<RooRealVar*> histVars(3, nullptr);
1219
1220 for(std::size_t iVar = 0; iVar < varNames.size(); ++iVar) {
1221 if(varNames[iVar].empty()) continue;
1222 if(iVar >= 3) {
1223 std::stringstream errMsg;
1224 errMsg << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR more than three variable names passed, but maximum number of supported variables is three";
1225 coutE(Plotting) << errMsg.str() << std::endl;
1226 throw std::invalid_argument(errMsg.str());
1227 }
1228 auto var = static_cast<RooRealVar*>(vars->find(varNames[iVar].c_str()));
1229 if(!var) {
1230 std::stringstream errMsg;
1231 errMsg << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR variable " << varNames[iVar] << " does not exist in argset: " << *vars;
1232 coutE(Plotting) << errMsg.str() << std::endl;
1233 throw std::runtime_error(errMsg.str());
1234 }
1235 histVars[iVar] = var;
1236 }
1237
1238 // Construct list of named arguments to pass to the implementation version of createHistogram()
1239
1240 RooLinkedList argList ;
1241 if (xbins>0) {
1242 argList.Add(RooFit::Binning(xbins).Clone()) ;
1243 }
1244
1245 if (histVars[1]) {
1246 argList.Add(RooFit::YVar(*histVars[1], ybins > 0 ? RooFit::Binning(ybins) : RooCmdArg::none()).Clone()) ;
1247 }
1248
1249 if (histVars[2]) {
1250 argList.Add(RooFit::ZVar(*histVars[2], zbins > 0 ? RooFit::Binning(zbins) : RooCmdArg::none()).Clone()) ;
1251 }
1252
1253 // Call implementation function
1254 TH1* result = createHistogram(GetName(), *histVars[0], argList) ;
1255
1256 // Delete temporary list of RooCmdArgs
1257 argList.Delete() ;
1258
1259 return result ;
1260}
1261
1262
1263
1264////////////////////////////////////////////////////////////////////////////////
1265/// Create and fill a ROOT histogram TH1, TH2 or TH3 with the values of this function.
1266///
1267/// \param[in] name Name of the ROOT histogram
1268/// \param[in] xvar Observable to be std::mapped on x axis of ROOT histogram
1269/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 Arguments according to list below
1270/// \return TH1 *, one of TH{1,2,3}. The caller takes ownership.
1271///
1272/// <table>
1273/// <tr><th><th> Effect on histogram creation
1274/// <tr><td> `IntrinsicBinning()` <td> Apply binning defined by function or pdf (as advertised via binBoundaries() method)
1275/// <tr><td> `Binning(const char* name)` <td> Apply binning with given name to x axis of histogram
1276/// <tr><td> `Binning(RooAbsBinning& binning)` <td> Apply specified binning to x axis of histogram
1277/// <tr><td> `Binning(int nbins, [double lo, double hi])` <td> Apply specified binning to x axis of histogram
1278/// <tr><td> `ConditionalObservables(Args_t &&... argsOrArgSet)` <td> Do not normalise PDF over following observables when projecting PDF into histogram.
1279// Arguments can either be multiple RooRealVar or a single RooArgSet containing them.
1280/// <tr><td> `Scaling(bool)` <td> Apply density-correction scaling (multiply by bin volume), default is true
1281/// <tr><td> `Extended(bool)` <td> Plot event yield instead of probability density (for extended pdfs only)
1282///
1283/// <tr><td> `YVar(const RooAbsRealLValue& var,...)` <td> Observable to be std::mapped on y axis of ROOT histogram.
1284/// The YVar() and ZVar() arguments can be supplied with optional Binning() arguments to control the binning of the Y and Z axes, e.g.
1285/// ```
1286/// createHistogram("histo",x,Binning(-1,1,20), YVar(y,Binning(-1,1,30)), ZVar(z,Binning("zbinning")))
1287/// ```
1288/// <tr><td> `ZVar(const RooAbsRealLValue& var,...)` <td> Observable to be std::mapped on z axis of ROOT histogram
1289/// </table>
1290///
1291///
1292
1294 const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
1295 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
1296{
1297
1299 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
1300 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
1301 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
1302 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
1303
1304 return createHistogram(name,xvar,l) ;
1305}
1306
1307
1308////////////////////////////////////////////////////////////////////////////////
1309/// Internal method implementing createHistogram
1310
1311TH1* RooAbsReal::createHistogram(const char *name, const RooAbsRealLValue& xvar, RooLinkedList& argList) const
1312{
1313
1314 // Define configuration for this method
1315 RooCmdConfig pc("RooAbsReal::createHistogram(" + std::string(GetName()) + ")");
1316 pc.defineInt("scaling","Scaling",0,1) ;
1317 pc.defineInt("intBinning","IntrinsicBinning",0,2) ;
1318 pc.defineInt("extended","Extended",0,2) ;
1319
1320 pc.defineSet("compSet","SelectCompSet",0);
1321 pc.defineString("compSpec","SelectCompSpec",0) ;
1322 pc.defineSet("projObs","ProjectedObservables",0,nullptr) ;
1323 pc.defineObject("yvar","YVar",0,nullptr) ;
1324 pc.defineObject("zvar","ZVar",0,nullptr) ;
1325 pc.defineMutex("SelectCompSet","SelectCompSpec") ;
1326 pc.defineMutex("IntrinsicBinning","Binning") ;
1327 pc.defineMutex("IntrinsicBinning","BinningName") ;
1328 pc.defineMutex("IntrinsicBinning","BinningSpec") ;
1329 pc.allowUndefined() ;
1330
1331 // Process & check varargs
1332 pc.process(argList) ;
1333 if (!pc.ok(true)) {
1334 return nullptr ;
1335 }
1336
1337 RooArgList vars(xvar) ;
1338 RooAbsArg* yvar = static_cast<RooAbsArg*>(pc.getObject("yvar")) ;
1339 if (yvar) {
1340 vars.add(*yvar) ;
1341 }
1342 RooAbsArg* zvar = static_cast<RooAbsArg*>(pc.getObject("zvar")) ;
1343 if (zvar) {
1344 vars.add(*zvar) ;
1345 }
1346
1347 auto projObs = pc.getSet("projObs");
1348 RooArgSet* intObs = nullptr ;
1349
1350 bool doScaling = pc.getInt("scaling") ;
1351 Int_t doIntBinning = pc.getInt("intBinning") ;
1352 Int_t doExtended = pc.getInt("extended") ;
1353
1354 // If doExtended is two, selection is automatic, set to 1 of pdf is extended, to zero otherwise
1355 const RooAbsPdf* pdfSelf = dynamic_cast<const RooAbsPdf*>(this) ;
1356 if (!pdfSelf && doExtended == 1) {
1357 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-pdf object, ignored" << std::endl ;
1358 doExtended=0 ;
1359 }
1360 if (pdfSelf && doExtended==1 && pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended) {
1361 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName() << ") WARNING extended mode requested for a non-extendable pdf, ignored" << std::endl ;
1362 doExtended=0 ;
1363 }
1364 if (pdfSelf && doExtended==2) {
1365 doExtended = pdfSelf->extendMode()==RooAbsPdf::CanNotBeExtended ? 0 : 1 ;
1366 } else if(!pdfSelf) {
1367 doExtended = 0;
1368 }
1369
1370 const char* compSpec = pc.getString("compSpec") ;
1371 const RooArgSet* compSet = pc.getSet("compSet");
1372 bool haveCompSel = ( (compSpec && strlen(compSpec)>0) || compSet) ;
1373
1374 std::unique_ptr<RooBinning> intBinning;
1375 if (doIntBinning>0) {
1376 // Given RooAbsPdf* pdf and RooRealVar* obs
1377 std::unique_ptr<std::list<double>> bl{binBoundaries(const_cast<RooAbsRealLValue&>(xvar),xvar.getMin(),xvar.getMax())};
1378 if (!bl) {
1379 // Only emit warning when intrinsic binning is explicitly requested
1380 if (doIntBinning==1) {
1381 coutW(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
1382 << ") WARNING, intrinsic model binning requested for histogram, but model does not define bin boundaries, reverting to default binning"<< std::endl ;
1383 }
1384 } else {
1385 if (doIntBinning==2) {
1386 coutI(InputArguments) << "RooAbsReal::createHistogram(" << GetName()
1387 << ") INFO: Model has intrinsic binning definition, selecting that binning for the histogram"<< std::endl ;
1388 }
1389 std::vector<double> edges(bl->size());
1390 int i=0 ;
1391 for (auto const& elem : *bl) { edges[i++] = elem ; }
1392 intBinning = std::make_unique<RooBinning>(bl->size()-1,edges.data()) ;
1393 }
1394 }
1395
1396 RooLinkedList argListCreate(argList) ;
1397 RooCmdConfig::stripCmdList(argListCreate,"Scaling,ProjectedObservables,IntrinsicBinning,SelectCompSet,SelectCompSpec,Extended") ;
1398
1399 TH1* histo(nullptr) ;
1400 if (intBinning) {
1401 RooCmdArg tmp = RooFit::Binning(*intBinning) ;
1402 argListCreate.Add(&tmp) ;
1403 histo = xvar.createHistogram(name,argListCreate) ;
1404 } else {
1405 histo = xvar.createHistogram(name,argListCreate) ;
1406 }
1407
1408 // Do component selection here
1409 if (haveCompSel) {
1410
1411 // Get complete set of tree branch nodes
1412 RooArgSet branchNodeSet ;
1413 branchNodeServerList(&branchNodeSet) ;
1414
1415 // Discard any non-RooAbsReal nodes
1416 for(RooAbsArg * arg : branchNodeSet) {
1417 if (!dynamic_cast<RooAbsReal*>(arg)) {
1418 branchNodeSet.remove(*arg) ;
1419 }
1420 }
1421
1422 std::unique_ptr<RooArgSet> dirSelNodes;
1423 if (compSet) {
1424 dirSelNodes.reset(static_cast<RooArgSet*>(branchNodeSet.selectCommon(*compSet)));
1425 } else {
1426 dirSelNodes.reset(static_cast<RooArgSet*>(branchNodeSet.selectByName(compSpec)));
1427 }
1428 if (!dirSelNodes->empty()) {
1429 coutI(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") directly selected PDF components: " << *dirSelNodes << std::endl ;
1430
1431 // Do indirect selection and activate both
1432 plotOnCompSelect(dirSelNodes.get()) ;
1433 } else {
1434 if (compSet) {
1435 coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection set " << *compSet << " does not match any components of p.d.f." << std::endl ;
1436 } else {
1437 coutE(Plotting) << "RooAbsPdf::createHistogram(" << GetName() << ") ERROR: component selection expression '" << compSpec << "' does not select any components of p.d.f." << std::endl ;
1438 }
1439 return nullptr ;
1440 }
1441 }
1442
1443 double scaleFactor(1.0) ;
1444 if (doExtended) {
1445 scaleFactor = pdfSelf->expectedEvents(vars) ;
1446 doScaling=false ;
1447 }
1448
1449 fillHistogram(histo,vars,scaleFactor,intObs,doScaling,projObs,false) ;
1450
1451 // Deactivate component selection
1452 if (haveCompSel) {
1453 plotOnCompSelect(nullptr) ;
1454 }
1455
1456
1457 return histo ;
1458}
1459
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Helper function for plotting of composite p.d.fs. Given
1463/// a set of selected components that should be plotted,
1464/// find all nodes that (in)directly depend on these selected
1465/// nodes. Mark all directly and indirectly selected nodes
1466/// as 'selected' using the selectComp() method
1467
1469{
1470 // Get complete set of tree branch nodes
1471 RooArgSet branchNodeSet;
1472 branchNodeServerList(&branchNodeSet);
1473
1474 // Discard any non-PDF nodes
1475 // Iterate by number because collection is being modified! Iterators may invalidate ...
1476 for (unsigned int i = 0; i < branchNodeSet.size(); ++i) {
1477 const auto arg = branchNodeSet[i];
1478 if (!dynamic_cast<RooAbsReal*>(arg)) {
1479 branchNodeSet.remove(*arg) ;
1480 }
1481 }
1482
1483 // If no set is specified, restored all selection bits to true
1484 if (!selNodes) {
1485 // Reset PDF selection bits to true
1486 for (const auto arg : branchNodeSet) {
1487 static_cast<RooAbsReal*>(arg)->selectComp(true);
1488 }
1489 return ;
1490 }
1491
1492
1493 // Add all nodes below selected nodes that are value servers
1494 RooArgSet tmp;
1495 for (const auto arg : branchNodeSet) {
1496 for (const auto selNode : *selNodes) {
1497 if (selNode->dependsOn(*arg, nullptr, /*valueOnly=*/true)) {
1498 tmp.add(*arg,true);
1499 }
1500 }
1501 }
1502
1503 // Add all nodes that depend on selected nodes by value
1504 for (const auto arg : branchNodeSet) {
1505 if (arg->dependsOn(*selNodes, nullptr, /*valueOnly=*/true)) {
1506 tmp.add(*arg,true);
1507 }
1508 }
1509
1510 tmp.remove(*selNodes, true);
1511 tmp.remove(*this);
1512 selNodes->add(tmp);
1513 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") indirectly selected PDF components: " << tmp << std::endl ;
1514
1515 // Set PDF selection bits according to selNodes
1516 for (const auto arg : branchNodeSet) {
1517 bool select = selNodes->find(arg->GetName()) != nullptr;
1518 static_cast<RooAbsReal*>(arg)->selectComp(select);
1519 }
1520}
1521
1522
1523
1524////////////////////////////////////////////////////////////////////////////////
1525/// Plot (project) PDF on specified frame. If a PDF is plotted in an empty frame, it
1526/// will show a unit normalized curve in the frame variable, taken at the present value
1527/// of other observables defined for this PDF.
1528///
1529/// \param[in] frame pointer to RooPlot
1530/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9,arg10 Ordered arguments
1531///
1532/// If a PDF is plotted in a frame in which a dataset has already been plotted, it will
1533/// show a projected curve integrated over all variables that were present in the shown
1534/// dataset except for the one on the x-axis. The normalization of the curve will also
1535/// be adjusted to the event count of the plotted dataset. An informational message
1536/// will be printed for each projection step that is performed.
1537///
1538/// This function takes the following named arguments
1539/// <table>
1540/// <tr><th><th> Projection control
1541/// <tr><td> `Slice(const RooArgSet& set)` <td> Override default projection behaviour by omitting observables listed
1542/// in set from the projection, i.e. by not integrating over these.
1543/// Slicing is usually only sensible in discrete observables, by e.g. creating a slice
1544/// of the PDF at the current value of the category observable.
1545///
1546/// <tr><td> `Slice(RooCategory& cat, const char* label)` <td> Override default projection behaviour by omitting the specified category
1547/// observable from the projection, i.e., by not integrating over all states of this category.
1548/// The slice is positioned at the given label value. To pass multiple Slice() commands, please use the
1549/// Slice(std::map<RooCategory*, std::string> const&) argument explained below.
1550///
1551/// <tr><td> `Slice(std::map<RooCategory*, std::string> const&)` <td> Omits multiple categories from the projection, as explianed above.
1552/// Can be used with initializer lists for convenience, e.g.
1553/// ```{.cpp}
1554/// pdf.plotOn(frame, Slice({{&tagCategory, "2tag"}, {&jetCategory, "3jet"}});
1555/// ```
1556///
1557/// <tr><td> `Project(const RooArgSet& set)` <td> Override default projection behaviour by projecting over observables
1558/// given in the set, ignoring the default projection behavior. Advanced use only.
1559///
1560/// <tr><td> `ProjWData(const RooAbsData& d)` <td> Override default projection _technique_ (integration). For observables present in given dataset
1561/// projection of PDF is achieved by constructing an average over all observable values in given set.
1562/// Consult RooFit plotting tutorial for further explanation of meaning & use of this technique
1563///
1564/// <tr><td> `ProjWData(const RooArgSet& s, const RooAbsData& d)` <td> As above but only consider subset 's' of observables in dataset 'd' for projection through data averaging
1565///
1566/// <tr><td> `ProjectionRange(const char* rn)` <td> Override default range of projection integrals to a different range specified by given range name.
1567/// This technique allows you to project a finite width slice in a real-valued observable
1568///
1569/// <tr><td> `NumCPU(Int_t ncpu)` <td> Number of CPUs to use simultaneously to calculate data-weighted projections (only in combination with ProjWData)
1570///
1571///
1572/// <tr><th><th> Misc content control
1573/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per curve. A negative
1574/// value suppress output completely, a zero value will only print the error count per p.d.f component,
1575/// a positive value is will print details of each error up to numErr messages per p.d.f component.
1576///
1577/// <tr><td> `EvalErrorValue(double value)` <td> Set curve points at which (pdf) evaluation errors occur to specified value. By default the
1578/// function value is plotted.
1579///
1580/// <tr><td> `Normalization(double scale, ScaleType code)` <td> Adjust normalization by given scale factor. Interpretation of number depends on code:
1581/// - Relative: relative adjustment factor for a normalized function,
1582/// - NumEvent: scale to match given number of events.
1583/// - Raw: relative adjustment factor for an un-normalized function.
1584///
1585/// <tr><td> `Name(const chat* name)` <td> Give curve specified name in frame. Useful if curve is to be referenced later
1586///
1587/// <tr><td> `Asymmetry(const RooCategory& c)` <td> Show the asymmetry of the PDF in given two-state category [F(+)-F(-)] / [F(+)+F(-)] rather than
1588/// the PDF projection. Category must have two states with indices -1 and +1 or three states with
1589/// indices -1,0 and +1.
1590///
1591/// <tr><td> `ShiftToZero(bool flag)` <td> Shift entire curve such that lowest visible point is at exactly zero. Mostly useful when plotting \f$ -\log(L) \f$ or \f$ \chi^2 \f$ distributions
1592///
1593/// <tr><td> `AddTo(const char* name, double_t wgtSelf, double_t wgtOther)` <td> Add constructed projection to already existing curve with given name and relative weight factors
1594/// <tr><td> `Components(const char* names)` <td> When plotting sums of PDFs, plot only the named components (*e.g.* only
1595/// the signal of a signal+background model).
1596/// <tr><td> `Components(const RooArgSet& compSet)` <td> As above, but pass a RooArgSet of the components themselves.
1597///
1598/// <tr><th><th> Plotting control
1599/// <tr><td> `DrawOption(const char* opt)` <td> Select ROOT draw option for resulting TGraph object. Currently supported options are "F" (fill), "L" (line), and "P" (points).
1600/// \note Option "P" will cause RooFit to plot (and treat) this pdf as if it were data! This is intended for plotting "corrected data"-type pdfs such as "data-minus-background" or unfolded datasets.
1601///
1602/// <tr><td> `LineStyle(Int_t style)` <td> Select line style by ROOT line style code, default is solid
1603///
1604/// <tr><td> `LineColor(Int_t color)` <td> Select line color by ROOT color code, default is blue
1605///
1606/// <tr><td> `LineWidth(Int_t width)` <td> Select line with in pixels, default is 3
1607///
1608/// <tr><td> `MarkerStyle(Int_t style)` <td> Select the ROOT marker style, default is 21
1609///
1610/// <tr><td> `MarkerColor(Int_t color)` <td> Select the ROOT marker color, default is black
1611///
1612/// <tr><td> `MarkerSize(double size)` <td> Select the ROOT marker size
1613///
1614/// <tr><td> `FillStyle(Int_t style)` <td> Select fill style, default is not filled. If a filled style is selected, also use VLines()
1615/// to add vertical downward lines at end of curve to ensure proper closure. Add `DrawOption("F")` for filled drawing.
1616/// <tr><td> `FillColor(Int_t color)` <td> Select fill color by ROOT color code
1617///
1618/// <tr><td> `Range(const char* name)` <td> Only draw curve in range defined by given name
1619///
1620/// <tr><td> `Range(double lo, double hi)` <td> Only draw curve in specified range
1621///
1622/// <tr><td> `VLines()` <td> Add vertical lines to y=0 at end points of curve
1623///
1624/// <tr><td> `Precision(double eps)` <td> Control precision of drawn curve w.r.t to scale of plot, default is 1e-3. Higher precision
1625/// will result in more and more densely spaced curve points
1626///
1627/// <tr><td> `Invisible(bool flag)` <td> Add curve to frame, but do not display. Useful in combination AddTo()
1628///
1629/// <tr><td> `VisualizeError(const RooFitResult& fitres, double Z=1, bool linearMethod=true)`
1630/// <td> Visualize the uncertainty on the parameters, as given in fitres, at 'Z' sigma'. The linear method is fast but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made. Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations, but may take (much) longer to calculate
1631///
1632/// <tr><td> `VisualizeError(const RooFitResult& fitres, const RooArgSet& param, double Z=1, bool linearMethod=true)`
1633/// <td> Visualize the uncertainty on the subset of parameters 'param', as given in fitres, at 'Z' sigma'
1634/// </table>
1635///
1636/// Details on error band visualization
1637/// -----------------------------------
1638/// *VisualizeError() uses plotOnWithErrorBand(). Documentation of the latter:*
1639/// \see plotOnWithErrorBand()
1640
1641RooPlot* RooAbsReal::plotOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2,
1642 const RooCmdArg& arg3, const RooCmdArg& arg4,
1643 const RooCmdArg& arg5, const RooCmdArg& arg6,
1644 const RooCmdArg& arg7, const RooCmdArg& arg8,
1645 const RooCmdArg& arg9, const RooCmdArg& arg10) const
1646{
1648 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
1649 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
1650 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
1651 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
1652 l.Add((TObject*)&arg9) ; l.Add((TObject*)&arg10) ;
1653 return plotOn(frame,l) ;
1654}
1655
1656
1657
1658////////////////////////////////////////////////////////////////////////////////
1659/// Internal back-end function of plotOn() with named arguments
1660
1662{
1663 // Special handling here if argList contains RangeWithName argument with multiple
1664 // range names -- Need to translate this call into multiple calls
1665
1666 RooCmdArg* rcmd = (RooCmdArg*) argList.FindObject("RangeWithName") ;
1667 if (rcmd && TString(rcmd->getString(0)).Contains(",")) {
1668
1669 // List joint ranges as choice of normalization for all later processing
1670 RooCmdArg rnorm = RooFit::NormRange(rcmd->getString(0)) ;
1671 argList.Add(&rnorm) ;
1672
1673 std::vector<std::string> rlist;
1674
1675 // Separate named ranges using strtok
1676 for (const std::string& rangeNameToken : ROOT::Split(rcmd->getString(0), ",")) {
1677 rlist.emplace_back(rangeNameToken);
1678 }
1679
1680 for (const auto& rangeString : rlist) {
1681 // Process each range with a separate command with a single range to be plotted
1682 rcmd->setString(0, rangeString.c_str());
1683 RooAbsReal::plotOn(frame,argList);
1684 }
1685 return frame ;
1686
1687 }
1688
1689 // Define configuration for this method
1690 RooCmdConfig pc("RooAbsReal::plotOn(" + std::string(GetName()) + ")");
1691 pc.defineString("drawOption","DrawOption",0,"L") ;
1692 pc.defineString("projectionRangeName","ProjectionRange",0,"",true) ;
1693 pc.defineString("curveNameSuffix","CurveNameSuffix",0,"") ;
1694 pc.defineString("sliceCatState","SliceCat",0,"",true) ;
1695 pc.defineDouble("scaleFactor","Normalization",0,1.0) ;
1696 pc.defineInt("scaleType","Normalization",0,Relative) ;
1697 pc.defineSet("sliceSet","SliceVars",0) ;
1698 pc.defineObject("sliceCatList","SliceCat",0,nullptr,true) ;
1699 // This dummy is needed for plotOn to recognize the "SliceCatMany" command.
1700 // It is not used directly, but the "SliceCat" commands are nested in it.
1701 // Removing this dummy definition results in "ERROR: unrecognized command: SliceCatMany".
1702 pc.defineObject("dummy1","SliceCatMany",0) ;
1703 pc.defineSet("projSet","Project",0) ;
1704 pc.defineObject("asymCat","Asymmetry",0) ;
1705 pc.defineDouble("precision","Precision",0,1e-3) ;
1706 pc.defineDouble("evalErrorVal","EvalErrorValue",0,0) ;
1707 pc.defineInt("doEvalError","EvalErrorValue",0,0) ;
1708 pc.defineInt("shiftToZero","ShiftToZero",0,0) ;
1709 pc.defineSet("projDataSet","ProjData",0) ;
1710 pc.defineObject("projData","ProjData",1) ;
1711 pc.defineObject("errorFR","VisualizeError",0) ;
1712 pc.defineDouble("errorZ","VisualizeError",0,1.) ;
1713 pc.defineSet("errorPars","VisualizeError",0) ;
1714 pc.defineInt("linearMethod","VisualizeError",0,0) ;
1715 pc.defineInt("binProjData","ProjData",0,0) ;
1716 pc.defineDouble("rangeLo","Range",0,-999.) ;
1717 pc.defineDouble("rangeHi","Range",1,-999.) ;
1718 pc.defineInt("numee","PrintEvalErrors",0,10) ;
1719 pc.defineInt("rangeAdjustNorm","Range",0,0) ;
1720 pc.defineInt("rangeWNAdjustNorm","RangeWithName",0,0) ;
1721 pc.defineInt("VLines","VLines",0,2) ; // 2==ExtendedWings
1722 pc.defineString("rangeName","RangeWithName",0,"") ;
1723 pc.defineString("normRangeName","NormRange",0,"") ;
1724 pc.defineInt("markerColor","MarkerColor",0,-999) ;
1725 pc.defineInt("markerStyle","MarkerStyle",0,-999) ;
1726 pc.defineDouble("markerSize","MarkerSize",0,-999) ;
1727 pc.defineInt("lineColor","LineColor",0,-999) ;
1728 pc.defineInt("lineStyle","LineStyle",0,-999) ;
1729 pc.defineInt("lineWidth","LineWidth",0,-999) ;
1730 pc.defineInt("fillColor","FillColor",0,-999) ;
1731 pc.defineInt("fillStyle","FillStyle",0,-999) ;
1732 pc.defineString("curveName","Name",0,"") ;
1733 pc.defineInt("curveInvisible","Invisible",0,0) ;
1734 pc.defineInt("showProg","ShowProgress",0,0) ;
1735 pc.defineInt("numCPU","NumCPU",0,1) ;
1736 pc.defineInt("interleave","NumCPU",1,0) ;
1737 pc.defineString("addToCurveName","AddTo",0,"") ;
1738 pc.defineDouble("addToWgtSelf","AddTo",0,1.) ;
1739 pc.defineDouble("addToWgtOther","AddTo",1,1.) ;
1740 pc.defineInt("moveToBack","MoveToBack",0,0) ;
1741 pc.defineMutex("SliceVars","Project") ;
1742 pc.defineMutex("AddTo","Asymmetry") ;
1743 pc.defineMutex("Range","RangeWithName") ;
1744 pc.defineMutex("VisualizeError","VisualizeErrorData") ;
1745
1746 // Process & check varargs
1747 pc.process(argList) ;
1748 if (!pc.ok(true)) {
1749 return frame ;
1750 }
1751
1752 PlotOpt o ;
1753 TString drawOpt(pc.getString("drawOption"));
1754
1755 RooFitResult* errFR = (RooFitResult*) pc.getObject("errorFR") ;
1756 double errZ = pc.getDouble("errorZ") ;
1757 RooArgSet* errPars = pc.getSet("errorPars") ;
1758 bool linMethod = pc.getInt("linearMethod") ;
1759 if (!drawOpt.Contains("P") && errFR) {
1760 return plotOnWithErrorBand(frame,*errFR,errZ,errPars,argList,linMethod) ;
1761 } else {
1762 o.errorFR = errFR;
1763 }
1764
1765 // Extract values from named arguments
1766 o.numee = pc.getInt("numee") ;
1767 o.drawOptions = drawOpt.Data();
1768 o.curveNameSuffix = pc.getString("curveNameSuffix") ;
1769 o.scaleFactor = pc.getDouble("scaleFactor") ;
1770 o.stype = (ScaleType) pc.getInt("scaleType") ;
1771 o.projData = (const RooAbsData*) pc.getObject("projData") ;
1772 o.binProjData = pc.getInt("binProjData") ;
1773 o.projDataSet = pc.getSet("projDataSet");
1774 o.numCPU = pc.getInt("numCPU") ;
1775 o.interleave = (RooFit::MPSplit) pc.getInt("interleave") ;
1776 o.eeval = pc.getDouble("evalErrorVal") ;
1777 o.doeeval = pc.getInt("doEvalError") ;
1778
1779 const RooArgSet* sliceSetTmp = pc.getSet("sliceSet");
1780 std::unique_ptr<RooArgSet> sliceSet{sliceSetTmp ? static_cast<RooArgSet*>(sliceSetTmp->Clone()) : nullptr};
1781 const RooArgSet* projSet = pc.getSet("projSet") ;
1782 const RooAbsCategoryLValue* asymCat = (const RooAbsCategoryLValue*) pc.getObject("asymCat") ;
1783
1784
1785 // Look for category slice arguments and add them to the master slice list if found
1786 const char* sliceCatState = pc.getString("sliceCatState",nullptr,true) ;
1787 const RooLinkedList& sliceCatList = pc.getObjectList("sliceCatList") ;
1788 if (sliceCatState) {
1789
1790 // Make the master slice set if it doesnt exist
1791 if (!sliceSet) {
1792 sliceSet = std::make_unique<RooArgSet>();
1793 }
1794
1795 // Loop over all categories provided by (multiple) Slice() arguments
1796 auto catTokens = ROOT::Split(sliceCatState, ",");
1797 auto iter = sliceCatList.begin();
1798 for (unsigned int i=0; i < catTokens.size(); ++i) {
1799 if (auto scat = static_cast<RooCategory*>(*iter)) {
1800 // Set the slice position to the value indicate by slabel
1801 scat->setLabel(catTokens[i]) ;
1802 // Add the slice category to the master slice set
1803 sliceSet->add(*scat,false) ;
1804 }
1805 ++iter;
1806 }
1807 }
1808
1809 o.precision = pc.getDouble("precision") ;
1810 o.shiftToZero = (pc.getInt("shiftToZero")!=0) ;
1811 Int_t vlines = pc.getInt("VLines");
1812 if (pc.hasProcessed("Range")) {
1813 o.rangeLo = pc.getDouble("rangeLo") ;
1814 o.rangeHi = pc.getDouble("rangeHi") ;
1815 o.postRangeFracScale = pc.getInt("rangeAdjustNorm") ;
1816 if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
1817 } else if (pc.hasProcessed("RangeWithName")) {
1818 o.normRangeName = pc.getString("rangeName",nullptr,true) ;
1819 o.rangeLo = frame->getPlotVar()->getMin(pc.getString("rangeName",nullptr,true)) ;
1820 o.rangeHi = frame->getPlotVar()->getMax(pc.getString("rangeName",nullptr,true)) ;
1821 o.postRangeFracScale = pc.getInt("rangeWNAdjustNorm") ;
1822 if (vlines==2) vlines=0 ; // Default is NoWings if range was specified
1823 }
1824
1825
1826 // If separate normalization range was specified this overrides previous settings
1827 if (pc.hasProcessed("NormRange")) {
1828 o.normRangeName = pc.getString("normRangeName") ;
1829 o.postRangeFracScale = true ;
1830 }
1831
1832 o.wmode = (vlines==2)?RooCurve::Extended:(vlines==1?RooCurve::Straight:RooCurve::NoWings) ;
1833 o.projectionRangeName = pc.getString("projectionRangeName",nullptr,true) ;
1834 o.curveName = pc.getString("curveName",nullptr,true) ;
1835 o.curveInvisible = pc.getInt("curveInvisible") ;
1836 o.progress = pc.getInt("showProg") ;
1837 o.addToCurveName = pc.getString("addToCurveName",nullptr,true) ;
1838 o.addToWgtSelf = pc.getDouble("addToWgtSelf") ;
1839 o.addToWgtOther = pc.getDouble("addToWgtOther") ;
1840
1842 coutE(InputArguments) << "RooAbsReal::plotOn(" << GetName() << ") cannot find existing curve " << o.addToCurveName << " to add to in RooPlot" << std::endl ;
1843 return frame ;
1844 }
1845
1846 RooArgSet projectedVars ;
1847 if (sliceSet) {
1848 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have slice " << *sliceSet << std::endl ;
1849
1850 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1851
1852 // Take out the sliced variables
1853 for (const auto sliceArg : *sliceSet) {
1854 RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
1855 if (arg) {
1856 projectedVars.remove(*arg) ;
1857 } else {
1858 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") slice variable "
1859 << sliceArg->GetName() << " was not projected anyway" << std::endl ;
1860 }
1861 }
1862 } else if (projSet) {
1863 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have projSet " << *projSet << std::endl ;
1864 makeProjectionSet(frame->getPlotVar(),projSet,projectedVars,false) ;
1865 } else {
1866 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: have neither sliceSet nor projSet " << std::endl ;
1867 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1868 }
1869 o.projSet = &projectedVars ;
1870
1871 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") Preprocessing: projectedVars = " << projectedVars << std::endl ;
1872
1873
1874 RooPlot* ret ;
1875 if (!asymCat) {
1876 // Forward to actual calculation
1877 ret = RooAbsReal::plotOn(frame,o) ;
1878 } else {
1879 // Forward to actual calculation
1880 ret = RooAbsReal::plotAsymOn(frame,*asymCat,o) ;
1881 }
1882
1883 // Optionally adjust line/fill attributes
1884 Int_t lineColor = pc.getInt("lineColor") ;
1885 Int_t lineStyle = pc.getInt("lineStyle") ;
1886 Int_t lineWidth = pc.getInt("lineWidth") ;
1887 Int_t markerColor = pc.getInt("markerColor") ;
1888 Int_t markerStyle = pc.getInt("markerStyle") ;
1889 Size_t markerSize = pc.getDouble("markerSize") ;
1890 Int_t fillColor = pc.getInt("fillColor") ;
1891 Int_t fillStyle = pc.getInt("fillStyle") ;
1892 if (lineColor!=-999) ret->getAttLine()->SetLineColor(lineColor) ;
1893 if (lineStyle!=-999) ret->getAttLine()->SetLineStyle(lineStyle) ;
1894 if (lineWidth!=-999) ret->getAttLine()->SetLineWidth(lineWidth) ;
1895 if (fillColor!=-999) ret->getAttFill()->SetFillColor(fillColor) ;
1896 if (fillStyle!=-999) ret->getAttFill()->SetFillStyle(fillStyle) ;
1897 if (markerColor!=-999) ret->getAttMarker()->SetMarkerColor(markerColor) ;
1898 if (markerStyle!=-999) ret->getAttMarker()->SetMarkerStyle(markerStyle) ;
1899 if (markerSize!=-999) ret->getAttMarker()->SetMarkerSize(markerSize) ;
1900
1901 if ((fillColor != -999 || fillStyle != -999) && !drawOpt.Contains("F")) {
1902 coutW(Plotting) << "Fill color or style was set for plotting \"" << GetName()
1903 << "\", but these only have an effect when 'DrawOption(\"F\")' for fill is used at the same time." << std::endl;
1904 }
1905
1906 // Move last inserted object to back to drawing stack if requested
1907 if (pc.getInt("moveToBack") && frame->numItems()>1) {
1908 frame->drawBefore(frame->getObject(0)->GetName(), frame->getCurve()->GetName());
1909 }
1910
1911 return ret ;
1912}
1913
1914
1915
1916/// Plotting engine function for internal use
1917///
1918/// Plot ourselves on given frame. If frame contains a histogram, all dimensions of the plotted
1919/// function that occur in the previously plotted dataset are projected via partial integration,
1920/// otherwise no projections are performed. Optionally, certain projections can be performed
1921/// by summing over the values present in a provided dataset ('projData'), to correctly
1922/// project out data dependents that are not properly described by the PDF (e.g. per-event errors).
1923///
1924/// The functions value can be multiplied with an optional scale factor. The interpretation
1925/// of the scale factor is unique for generic real functions, for PDFs there are various interpretations
1926/// possible, which can be selection with 'stype' (see RooAbsPdf::plotOn() for details).
1927///
1928/// The default projection behaviour can be overridden by supplying an optional set of dependents
1929/// to project. For most cases, plotSliceOn() and plotProjOn() provide a more intuitive interface
1930/// to modify the default projection behaviour.
1931//_____________________________________________________________________________
1932// coverity[PASS_BY_VALUE]
1934{
1935
1936
1937 // Sanity checks
1938 if (plotSanityChecks(frame)) return frame ;
1939
1940 // ProjDataVars is either all projData observables, or the user indicated subset of it
1941 RooArgSet projDataVars ;
1942 if (o.projData) {
1943 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjData with observables = " << *o.projData->get() << std::endl ;
1944 if (o.projDataSet) {
1945 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(o.projData->get()->selectCommon(*o.projDataSet))};
1946 projDataVars.add(*tmp) ;
1947 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have ProjDataSet = " << *o.projDataSet << " will only use this subset of projData" << std::endl ;
1948 } else {
1949 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") using full ProjData" << std::endl ;
1950 projDataVars.add(*o.projData->get()) ;
1951 }
1952 }
1953
1954 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") ProjDataVars = " << projDataVars << std::endl ;
1955
1956 // Make list of variables to be projected
1957 RooArgSet projectedVars ;
1958 RooArgSet sliceSet ;
1959 if (o.projSet) {
1960 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") have input projSet = " << *o.projSet << std::endl ;
1961 makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,false) ;
1962 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") calculated projectedVars = " << *o.projSet << std::endl ;
1963
1964 // Print list of non-projected variables
1965 if (frame->getNormVars()) {
1966 RooArgSet sliceSetTmp;
1967 getObservables(frame->getNormVars(), sliceSetTmp) ;
1968
1969 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") frame->getNormVars() that are also observables = " << sliceSetTmp << std::endl ;
1970
1971 sliceSetTmp.remove(projectedVars,true,true) ;
1972 sliceSetTmp.remove(*frame->getPlotVar(),true,true) ;
1973
1974 if (o.projData) {
1975 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(projDataVars.selectCommon(*o.projSet))};
1976 sliceSetTmp.remove(*tmp,true,true) ;
1977 }
1978
1979 if (!sliceSetTmp.empty()) {
1980 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on "
1981 << frame->getPlotVar()->GetName() << " represents a slice in " << sliceSetTmp << std::endl ;
1982 }
1983 sliceSet.add(sliceSetTmp) ;
1984 }
1985 } else {
1986 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
1987 }
1988
1989 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") projectedVars = " << projectedVars << " sliceSet = " << sliceSet << std::endl ;
1990
1991
1992 RooArgSet* projDataNeededVars = nullptr ;
1993 // Take out data-projected dependents from projectedVars
1994 if (o.projData) {
1995 projDataNeededVars = (RooArgSet*) projectedVars.selectCommon(projDataVars) ;
1996 projectedVars.remove(projDataVars,true,true) ;
1997 }
1998
1999 // Get the plot variable and remember its original value
2000 auto* plotVar = static_cast<RooRealVar*>(frame->getPlotVar());
2001 double oldPlotVarVal = plotVar->getVal();
2002
2003 // Inform user about projections
2004 if (!projectedVars.empty()) {
2005 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2006 << " integrates over variables " << projectedVars
2007 << (o.projectionRangeName?Form(" in range %s",o.projectionRangeName):"") << std::endl;
2008 }
2009 if (projDataNeededVars && !projDataNeededVars->empty()) {
2010 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2011 << " averages using data variables " << *projDataNeededVars << std::endl ;
2012 }
2013
2014 // Create projection integral
2015 RooArgSet* projectionCompList = nullptr ;
2016
2017 RooArgSet deps;
2018 getObservables(frame->getNormVars(), deps) ;
2019 deps.remove(projectedVars,true,true) ;
2020 if (projDataNeededVars) {
2021 deps.remove(*projDataNeededVars,true,true) ;
2022 }
2023 deps.remove(*plotVar,true,true) ;
2024 deps.add(*plotVar) ;
2025
2026 // Now that we have the final set of dependents, call checkObservables()
2027
2028 // WVE take out conditional observables
2029 if (checkObservables(&deps)) {
2030 coutE(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") error in checkObservables, abort" << std::endl ;
2031 if (projDataNeededVars) delete projDataNeededVars ;
2032 return frame ;
2033 }
2034
2035 RooAbsReal *projection = (RooAbsReal*) createPlotProjection(deps, &projectedVars, projectionCompList, o.projectionRangeName) ;
2036 cxcoutD(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot projection object is " << projection->GetName() << std::endl ;
2037 if (dologD(Plotting)) {
2038 projection->printStream(ccoutD(Plotting),0,kVerbose) ;
2039 }
2040
2041 // Always fix RooAddPdf normalizations
2042 RooArgSet fullNormSet(deps) ;
2043 fullNormSet.add(projectedVars) ;
2044 if (projDataNeededVars && !projDataNeededVars->empty()) {
2045 fullNormSet.add(*projDataNeededVars) ;
2046 }
2047
2048 std::unique_ptr<RooArgSet> projectionComponents(projection->getComponents());
2049 for(auto * pdf : dynamic_range_cast<RooAbsPdf*>(*projectionComponents)) {
2050 if (pdf) {
2051 pdf->selectNormalization(&fullNormSet) ;
2052 }
2053 }
2054
2055 // Apply data projection, if requested
2056 if (o.projData && projDataNeededVars && !projDataNeededVars->empty()) {
2057
2058 // If data set contains more rows than needed, make reduced copy first
2059 RooAbsData* projDataSel = (RooAbsData*)o.projData;
2060 std::unique_ptr<RooAbsData> projDataSelOwned;
2061
2062 if (projDataNeededVars->size() < o.projData->get()->size()) {
2063
2064 // Determine if there are any slice variables in the projection set
2065 std::unique_ptr<RooArgSet> sliceDataSet{static_cast<RooArgSet*>(sliceSet.selectCommon(*o.projData->get()))};
2066 TString cutString ;
2067 if (!sliceDataSet->empty()) {
2068 bool first(true) ;
2069 for(RooAbsArg * sliceVar : *sliceDataSet) {
2070 if (!first) {
2071 cutString.Append("&&") ;
2072 } else {
2073 first=false ;
2074 }
2075
2076 RooAbsRealLValue* real ;
2078 if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
2079 cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
2080 } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
2081 cutString.Append(Form("%s==%d",cat->GetName(),cat->getCurrentIndex())) ;
2082 }
2083 }
2084 }
2085
2086 if (!cutString.IsNull()) {
2087 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") reducing given projection dataset to entries with " << cutString << std::endl ;
2088 }
2089 projDataSelOwned = std::unique_ptr<RooAbsData>{const_cast<RooAbsData*>(o.projData)->reduce(*projDataNeededVars, cutString.IsNull() ? nullptr : cutString)};
2090 projDataSel = projDataSelOwned.get();
2091 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName()
2092 << ") only the following components of the projection data will be used: " << *projDataNeededVars << std::endl ;
2093 }
2094
2095 // Request binning of unbinned projection dataset that consists exclusively of category observables
2096 if (!o.binProjData && dynamic_cast<RooDataSet*>(projDataSel)!=nullptr) {
2097
2098 // Determine if dataset contains only categories
2099 bool allCat(true) ;
2100 for(RooAbsArg * arg2 : *projDataSel->get()) {
2101 if (!dynamic_cast<RooCategory*>(arg2)) allCat = false ;
2102 }
2103 if (allCat) {
2104 o.binProjData = true ;
2105 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") unbinned projection dataset consist only of discrete variables,"
2106 << " performing projection with binned copy for optimization." << std::endl ;
2107
2108 }
2109 }
2110
2111 // Bin projection dataset if requested
2112 if (o.binProjData) {
2113 projDataSelOwned = std::make_unique<RooDataHist>(std::string(projDataSel->GetName()) + "_binned","Binned projection data",*projDataSel->get(),*projDataSel);
2114 projDataSel = projDataSelOwned.get();
2115 }
2116
2117 // Construct scaled data weighted average
2118 ScaledDataWeightedAverage scaleBind{*projection, *projDataSel, o.scaleFactor, *plotVar};
2119
2120 // Set default range, if not specified
2121 if (o.rangeLo==0 && o.rangeHi==0) {
2122 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2123 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2124 }
2125
2126 // Construct name of curve for data weighed average
2127 std::string curveName(projection->GetName()) ;
2128 curveName.append("_DataAvg[" + projDataSel->get()->contentsString() + "]");
2129 // Append slice set specification if any
2130 if (!sliceSet.empty()) {
2131 curveName.append("_Slice[" + sliceSet.contentsString() + "]");
2132 }
2133 // Append any suffixes imported from RooAbsPdf::plotOn
2134 if (o.curveNameSuffix) {
2135 curveName.append(o.curveNameSuffix) ;
2136 }
2137
2138 // Curve constructor for data weighted average
2140 RooCurve *curve = new RooCurve(projection->GetName(),projection->GetTitle(),scaleBind,
2143
2144 curve->SetName(curveName.c_str()) ;
2145
2146 // Add self to other curve if requested
2147 if (o.addToCurveName) {
2148 RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;
2149
2150 // Curve constructor for sum of curves
2151 RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
2152 sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
2153 delete curve ;
2154 curve = sumCurve ;
2155
2156 }
2157
2158 if (o.curveName) {
2159 curve->SetName(o.curveName) ;
2160 }
2161
2162 // add this new curve to the specified plot frame
2163 frame->addPlotable(curve, o.drawOptions, o.curveInvisible);
2164
2165 } else {
2166
2167 // Set default range, if not specified
2168 if (o.rangeLo==0 && o.rangeHi==0) {
2169 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2170 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2171 }
2172
2173 // Calculate a posteriori range fraction scaling if requested (2nd part of normalization correction for
2174 // result fit on subrange of data)
2175 if (o.postRangeFracScale) {
2176 if (!o.normRangeName) {
2177 o.normRangeName = "plotRange" ;
2178 plotVar->setRange("plotRange",o.rangeLo,o.rangeHi) ;
2179 }
2180
2181 // Evaluate fractional correction integral always on full p.d.f, not component.
2182 GlobalSelectComponentRAII selectCompRAII(true);
2183 std::unique_ptr<RooAbsReal> intFrac{projection->createIntegral(*plotVar,*plotVar,o.normRangeName)};
2184 if(o.stype != RooAbsReal::Raw || this->InheritsFrom(RooAbsPdf::Class())){
2185 // this scaling should only be !=1 when plotting partial ranges
2186 // still, raw means raw
2187 o.scaleFactor /= intFrac->getVal() ;
2188 }
2189 }
2190
2191 // create a new curve of our function using the clone to do the evaluations
2192 // Curve constructor for regular projections
2193
2194 // Set default name of curve
2195 std::string curveName(projection->GetName()) ;
2196 if (!sliceSet.empty()) {
2197 curveName.append("_Slice[" + sliceSet.contentsString() + "]");
2198 }
2199 if (o.curveNameSuffix) {
2200 // Append any suffixes imported from RooAbsPdf::plotOn
2201 curveName.append(o.curveNameSuffix) ;
2202 }
2203
2204 TString opt(o.drawOptions);
2205 if(opt.Contains("P")){
2207 RooHist *graph= new RooHist(*projection,*plotVar,1.,o.scaleFactor,frame->getNormVars(),o.errorFR);
2209
2210 // Override name of curve by user name, if specified
2211 if (o.curveName) {
2212 graph->SetName(o.curveName) ;
2213 }
2214
2215 // add this new curve to the specified plot frame
2217 } else {
2219 RooCurve *curve = new RooCurve(*projection,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
2222 curve->SetName(curveName.c_str()) ;
2223
2224 // Add self to other curve if requested
2225 if (o.addToCurveName) {
2226 RooCurve* otherCurve = static_cast<RooCurve*>(frame->findObject(o.addToCurveName,RooCurve::Class())) ;
2227 RooCurve* sumCurve = new RooCurve(projection->GetName(),projection->GetTitle(),*curve,*otherCurve,o.addToWgtSelf,o.addToWgtOther) ;
2228 sumCurve->SetName(Form("%s_PLUS_%s",curve->GetName(),otherCurve->GetName())) ;
2229 delete curve ;
2230 curve = sumCurve ;
2231 }
2232
2233 // Override name of curve by user name, if specified
2234 if (o.curveName) {
2235 curve->SetName(o.curveName) ;
2236 }
2237
2238 // add this new curve to the specified plot frame
2239 frame->addPlotable(curve, o.drawOptions, o.curveInvisible);
2240 }
2241 }
2242
2243 if (projDataNeededVars) delete projDataNeededVars ;
2244 delete projectionCompList ;
2245 plotVar->setVal(oldPlotVarVal); // reset the plot variable value to not disturb the original state
2246 return frame;
2247}
2248
2249
2250
2251
2252////////////////////////////////////////////////////////////////////////////////
2253/// \deprecated OBSOLETE -- RETAINED FOR BACKWARD COMPATIBILITY. Use plotOn() with Slice() instead
2254
2255RooPlot* RooAbsReal::plotSliceOn(RooPlot *frame, const RooArgSet& sliceSet, Option_t* drawOptions,
2256 double scaleFactor, ScaleType stype, const RooAbsData* projData) const
2257{
2258 RooArgSet projectedVars ;
2259 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
2260
2261 // Take out the sliced variables
2262 for(RooAbsArg * sliceArg : sliceSet) {
2263 RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
2264 if (arg) {
2265 projectedVars.remove(*arg) ;
2266 } else {
2267 coutI(Plotting) << "RooAbsReal::plotSliceOn(" << GetName() << ") slice variable "
2268 << sliceArg->GetName() << " was not projected anyway" << std::endl ;
2269 }
2270 }
2271
2272 PlotOpt o ;
2273 o.drawOptions = drawOptions ;
2274 o.scaleFactor = scaleFactor ;
2275 o.stype = stype ;
2276 o.projData = projData ;
2277 o.projSet = &projectedVars ;
2278 return plotOn(frame,o) ;
2279}
2280
2281
2282
2283
2284//_____________________________________________________________________________
2285// coverity[PASS_BY_VALUE]
2287
2288{
2289 // Plotting engine for asymmetries. Implements the functionality if plotOn(frame,Asymmetry(...)))
2290 //
2291 // Plot asymmetry of ourselves, defined as
2292 //
2293 // asym = f(asymCat=-1) - f(asymCat=+1) / ( f(asymCat=-1) + f(asymCat=+1) )
2294 //
2295 // on frame. If frame contains a histogram, all dimensions of the plotted
2296 // asymmetry function that occur in the previously plotted dataset are projected via partial integration.
2297 // Otherwise no projections are performed,
2298 //
2299 // The asymmetry function can be multiplied with an optional scale factor. The default projection
2300 // behaviour can be overridden by supplying an optional set of dependents to project.
2301
2302 // Sanity checks
2303 if (plotSanityChecks(frame)) return frame ;
2304
2305 // ProjDataVars is either all projData observables, or the user indicated subset of it
2306 RooArgSet projDataVars ;
2307 if (o.projData) {
2308 if (o.projDataSet) {
2309 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(o.projData->get()->selectCommon(*o.projDataSet))};
2310 projDataVars.add(*tmp) ;
2311 } else {
2312 projDataVars.add(*o.projData->get()) ;
2313 }
2314 }
2315
2316 // Must depend on asymCat
2317 if (!dependsOn(asymCat)) {
2318 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2319 << ") function doesn't depend on asymmetry category " << asymCat.GetName() << std::endl ;
2320 return frame ;
2321 }
2322
2323 // asymCat must be a signCat
2324 if (!asymCat.isSignType()) {
2325 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2326 << ") asymmetry category must have 2 or 3 states with index values -1,0,1" << std::endl ;
2327 return frame ;
2328 }
2329
2330 // Make list of variables to be projected
2331 RooArgSet projectedVars ;
2332 RooArgSet sliceSet ;
2333 if (o.projSet) {
2334 makeProjectionSet(frame->getPlotVar(),o.projSet,projectedVars,false) ;
2335
2336 // Print list of non-projected variables
2337 if (frame->getNormVars()) {
2338 RooArgSet sliceSetTmp;
2339 getObservables(frame->getNormVars(), sliceSetTmp) ;
2340 sliceSetTmp.remove(projectedVars,true,true) ;
2341 sliceSetTmp.remove(*frame->getPlotVar(),true,true) ;
2342
2343 if (o.projData) {
2344 std::unique_ptr<RooArgSet> tmp{static_cast<RooArgSet*>(projDataVars.selectCommon(*o.projSet))};
2345 sliceSetTmp.remove(*tmp,true,true) ;
2346 }
2347
2348 if (!sliceSetTmp.empty()) {
2349 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on "
2350 << frame->getPlotVar()->GetName() << " represents a slice in " << sliceSetTmp << std::endl ;
2351 }
2352 sliceSet.add(sliceSetTmp) ;
2353 }
2354 } else {
2355 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
2356 }
2357
2358
2359 // Take out data-projected dependens from projectedVars
2360 RooArgSet* projDataNeededVars = nullptr ;
2361 if (o.projData) {
2362 projDataNeededVars = (RooArgSet*) projectedVars.selectCommon(projDataVars) ;
2363 projectedVars.remove(projDataVars,true,true) ;
2364 }
2365
2366 // Take out plotted asymmetry from projection
2367 if (projectedVars.find(asymCat.GetName())) {
2368 projectedVars.remove(*projectedVars.find(asymCat.GetName())) ;
2369 }
2370
2371 // Clone the plot variable
2372 RooAbsReal* realVar = (RooRealVar*) frame->getPlotVar() ;
2373 RooRealVar* plotVar = (RooRealVar*) realVar->Clone() ;
2374
2375 // Inform user about projections
2376 if (!projectedVars.empty()) {
2377 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") plot on " << plotVar->GetName()
2378 << " projects variables " << projectedVars << std::endl ;
2379 }
2380 if (projDataNeededVars && !projDataNeededVars->empty()) {
2381 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") plot on " << plotVar->GetName()
2382 << " averages using data variables "<< *projDataNeededVars << std::endl ;
2383 }
2384
2385
2386 // Customize two copies of projection with fixed negative and positive asymmetry
2387 std::unique_ptr<RooAbsCategoryLValue> asymPos{static_cast<RooAbsCategoryLValue*>(asymCat.Clone("asym_pos"))};
2388 std::unique_ptr<RooAbsCategoryLValue> asymNeg{static_cast<RooAbsCategoryLValue*>(asymCat.Clone("asym_neg"))};
2389 asymPos->setIndex(1) ;
2390 asymNeg->setIndex(-1) ;
2391 RooCustomizer custPos{*this,"pos"};
2392 RooCustomizer custNeg{*this,"neg"};
2393 //custPos->setOwning(true) ;
2394 //custNeg->setOwning(true) ;
2395 custPos.replaceArg(asymCat,*asymPos) ;
2396 custNeg.replaceArg(asymCat,*asymNeg) ;
2397 std::unique_ptr<RooAbsReal> funcPos{static_cast<RooAbsReal*>(custPos.build())};
2398 std::unique_ptr<RooAbsReal> funcNeg{static_cast<RooAbsReal*>(custNeg.build())};
2399
2400 // Create projection integral
2401 RooArgSet *posProjCompList, *negProjCompList ;
2402
2403 // Add projDataVars to normalized dependents of projection
2404 // This is needed only for asymmetries (why?)
2405 RooArgSet depPos(*plotVar,*asymPos) ;
2406 RooArgSet depNeg(*plotVar,*asymNeg) ;
2407 depPos.add(projDataVars) ;
2408 depNeg.add(projDataVars) ;
2409
2410 const RooAbsReal *posProj = funcPos->createPlotProjection(depPos, &projectedVars, posProjCompList, o.projectionRangeName) ;
2411 const RooAbsReal *negProj = funcNeg->createPlotProjection(depNeg, &projectedVars, negProjCompList, o.projectionRangeName) ;
2412 if (!posProj || !negProj) {
2413 coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() << ") Unable to create projections, abort" << std::endl ;
2414 return frame ;
2415 }
2416
2417 // Create a RooFormulaVar representing the asymmetry
2418 TString asymName(GetName()) ;
2419 asymName.Append("_Asym[") ;
2420 asymName.Append(asymCat.GetName()) ;
2421 asymName.Append("]") ;
2422 TString asymTitle(asymCat.GetName()) ;
2423 asymTitle.Append(" Asymmetry of ") ;
2424 asymTitle.Append(GetTitle()) ;
2425 RooFormulaVar funcAsym{asymName,asymTitle,"(@0-@1)/(@0+@1)",RooArgSet(*posProj,*negProj)};
2426
2427 if (o.projData) {
2428
2429 // If data set contains more rows than needed, make reduced copy first
2430 RooAbsData* projDataSel = (RooAbsData*)o.projData;
2431 std::unique_ptr<RooAbsData> projDataSelOwned;
2432 if (projDataNeededVars && projDataNeededVars->size() < o.projData->get()->size()) {
2433
2434 // Determine if there are any slice variables in the projection set
2435 RooArgSet sliceDataSet;
2436 sliceSet.selectCommon(*o.projData->get(), sliceDataSet);
2437 TString cutString ;
2438 if (!sliceDataSet.empty()) {
2439 bool first(true) ;
2440 for(RooAbsArg * sliceVar : sliceDataSet) {
2441 if (!first) {
2442 cutString.Append("&&") ;
2443 } else {
2444 first=false ;
2445 }
2446
2447 RooAbsRealLValue* real ;
2449 if ((real = dynamic_cast<RooAbsRealLValue*>(sliceVar))) {
2450 cutString.Append(Form("%s==%f",real->GetName(),real->getVal())) ;
2451 } else if ((cat = dynamic_cast<RooAbsCategoryLValue*>(sliceVar))) {
2452 cutString.Append(Form("%s==%d",cat->GetName(),cat->getCurrentIndex())) ;
2453 }
2454 }
2455 }
2456
2457 if (!cutString.IsNull()) {
2458 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2459 << ") reducing given projection dataset to entries with " << cutString << std::endl ;
2460 }
2461 projDataSelOwned = std::unique_ptr<RooAbsData>{const_cast<RooAbsData*>(o.projData)->reduce(*projDataNeededVars,cutString.IsNull() ? nullptr : cutString)};
2462 projDataSel = projDataSelOwned.get();
2463 coutI(Plotting) << "RooAbsReal::plotAsymOn(" << GetName()
2464 << ") only the following components of the projection data will be used: " << *projDataNeededVars << std::endl ;
2465 }
2466
2467
2468 // Construct scaled data weighted average
2469 ScaledDataWeightedAverage scaleBind{funcAsym, *projDataSel, o.scaleFactor, *plotVar};
2470
2471 // Set default range, if not specified
2472 if (o.rangeLo==0 && o.rangeHi==0) {
2473 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2474 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2475 }
2476
2477 // Construct name of curve for data weighed average
2478 TString curveName(funcAsym.GetName()) ;
2479 curveName.Append(Form("_DataAvg[%s]",projDataSel->get()->contentsString().c_str())) ;
2480 // Append slice set specification if any
2481 if (!sliceSet.empty()) {
2482 curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
2483 }
2484 // Append any suffixes imported from RooAbsPdf::plotOn
2485 if (o.curveNameSuffix) {
2486 curveName.Append(o.curveNameSuffix) ;
2487 }
2488
2489
2491 RooCurve *curve = new RooCurve(funcAsym.GetName(),funcAsym.GetTitle(),scaleBind,
2492 o.rangeLo,o.rangeHi,frame->GetNbinsX(),o.precision,o.precision,false,o.wmode,o.numee,o.doeeval,o.eeval) ;
2494
2495 dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;
2496 // add this new curve to the specified plot frame
2497 frame->addPlotable(curve, o.drawOptions);
2498
2499 ccoutW(Eval) << std::endl ;
2500 } else {
2501
2502 // Set default range, if not specified
2503 if (o.rangeLo==0 && o.rangeHi==0) {
2504 o.rangeLo = frame->GetXaxis()->GetXmin() ;
2505 o.rangeHi = frame->GetXaxis()->GetXmax() ;
2506 }
2507
2509 RooCurve* curve= new RooCurve(funcAsym,*plotVar,o.rangeLo,o.rangeHi,frame->GetNbinsX(),
2510 o.scaleFactor,nullptr,o.precision,o.precision,false,o.wmode,o.numee,o.doeeval,o.eeval);
2512
2513 dynamic_cast<TAttLine*>(curve)->SetLineColor(2) ;
2514
2515
2516 // Set default name of curve
2517 TString curveName(funcAsym.GetName()) ;
2518 if (!sliceSet.empty()) {
2519 curveName.Append(Form("_Slice[%s]",sliceSet.contentsString().c_str())) ;
2520 }
2521 if (o.curveNameSuffix) {
2522 // Append any suffixes imported from RooAbsPdf::plotOn
2523 curveName.Append(o.curveNameSuffix) ;
2524 }
2525 curve->SetName(curveName.Data()) ;
2526
2527 // add this new curve to the specified plot frame
2528 frame->addPlotable(curve, o.drawOptions);
2529
2530 }
2531
2532 // Cleanup
2533 delete posProjCompList ;
2534 delete negProjCompList ;
2535
2536 delete plotVar ;
2537
2538 return frame;
2539}
2540
2541
2542
2543////////////////////////////////////////////////////////////////////////////////
2544/// \brief Propagates parameter uncertainties to an uncertainty estimate for this RooAbsReal.
2545///
2546/// Estimates the uncertainty \f$\sigma_f(x;\theta)\f$ on a function \f$f(x;\theta)\f$ represented by this RooAbsReal.
2547/// Here, \f$\theta\f$ is a vector of parameters with uncertainties \f$\sigma_\theta\f$, and \f$x\f$ are usually observables.
2548/// The uncertainty is estimated by *linearly* propagating the parameter uncertainties using the correlation matrix from a fit result.
2549///
2550/// The square of the uncertainty on \f$f(x;\theta)\f$ is calculated as follows:
2551/// \f[
2552/// \sigma_f(x)^2 = \Delta f_i(x) \cdot \mathrm{Corr}_{i, j} \cdot \Delta f_j(x),
2553/// \f]
2554/// 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$
2555/// is the vector of function variations when changing the parameters one at a time, and
2556/// \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.
2557
2558double RooAbsReal::getPropagatedError(const RooFitResult &fr, const RooArgSet &nset) const
2559{
2560 // Calling getParameters() might be costly, but necessary to get the right
2561 // parameters in the RooAbsReal. The RooFitResult only stores snapshots.
2562 RooArgSet allParamsInAbsReal;
2563 getParameters(&nset, allParamsInAbsReal);
2564
2565 RooArgList paramList;
2566 for(auto * rrvFitRes : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {
2567
2568 auto rrvInAbsReal = static_cast<RooRealVar const*>(allParamsInAbsReal.find(*rrvFitRes));
2569
2570 // If this RooAbsReal is a RooRealVar in the fit result, we don't need to
2571 // propagate anything and can just return the error in the fit result
2572 if(rrvFitRes->namePtr() == namePtr()) return rrvFitRes->getError();
2573
2574 // Strip out parameters with zero error
2575 if (rrvFitRes->getError() <= std::abs(rrvFitRes->getVal()) * std::numeric_limits<double>::epsilon()) continue;
2576
2577 // Ignore parameters in the fit result that this RooAbsReal doesn't depend on
2578 if(!rrvInAbsReal) continue;
2579
2580 // Checking for float equality is a bad. We check if the values are
2581 // negligibly far away from each other, relative to the uncertainty.
2582 if(std::abs(rrvInAbsReal->getVal() - rrvFitRes->getVal()) > 0.01 * rrvFitRes->getError()) {
2583 std::stringstream errMsg;
2584 errMsg << "RooAbsReal::getPropagatedError(): the parameters of the RooAbsReal don't have"
2585 << " the same values as in the fit result! The logic of getPropagatedError is broken in this case.";
2586
2587 throw std::runtime_error(errMsg.str());
2588 }
2589
2590 paramList.add(*rrvInAbsReal);
2591 }
2592
2593 std::vector<double> plusVar;
2594 std::vector<double> minusVar;
2595 plusVar.reserve(paramList.size());
2596 minusVar.reserve(paramList.size());
2597
2598 // Create std::vector of plus,minus variations for each parameter
2599 TMatrixDSym V(paramList.size() == fr.floatParsFinal().size() ?
2600 fr.covarianceMatrix() :
2601 fr.reducedCovarianceMatrix(paramList)) ;
2602
2603 for (Int_t ivar=0 ; ivar<paramList.getSize() ; ivar++) {
2604
2605 auto& rrv = static_cast<RooRealVar&>(paramList[ivar]);
2606
2607 double cenVal = rrv.getVal() ;
2608 double errVal = sqrt(V(ivar,ivar)) ;
2609
2610 // Make Plus variation
2611 rrv.setVal(cenVal+errVal) ;
2612 plusVar.push_back(getVal(nset)) ;
2613
2614 // Make Minus variation
2615 rrv.setVal(cenVal-errVal) ;
2616 minusVar.push_back(getVal(nset)) ;
2617
2618 rrv.setVal(cenVal) ;
2619 }
2620
2621 // Re-evaluate this RooAbsReal with the central parameters just to be
2622 // extra-safe that a call to `getPropagatedError()` doesn't change any state.
2623 // It should not be necessary because thanks to the dirty flag propagation
2624 // the RooAbsReal is re-evaluated anyway the next time getVal() is called.
2625 // Still there are imaginable corner cases where it would not be triggered,
2626 // for example if the user changes the RooFit operation more after the error
2627 // propagation.
2628 getVal(nset);
2629
2630 TMatrixDSym C(paramList.getSize()) ;
2631 std::vector<double> errVec(paramList.size()) ;
2632 for (int i=0 ; i<paramList.getSize() ; i++) {
2633 errVec[i] = std::sqrt(V(i,i)) ;
2634 for (int j=i ; j<paramList.getSize() ; j++) {
2635 C(i,j) = V(i,j) / std::sqrt(V(i,i)*V(j,j));
2636 C(j,i) = C(i,j) ;
2637 }
2638 }
2639
2640 // Make std::vector of variations
2641 TVectorD F(plusVar.size()) ;
2642 for (unsigned int j=0 ; j<plusVar.size() ; j++) {
2643 F[j] = (plusVar[j]-minusVar[j])/2 ;
2644 }
2645
2646 // Calculate error in linear approximation from variations and correlation coefficient
2647 double sum = F*(C*F) ;
2648
2649 return sqrt(sum) ;
2650}
2651
2652
2653
2654////////////////////////////////////////////////////////////////////////////////
2655/// Plot function or PDF on frame with support for visualization of the uncertainty encoded in the given fit result fr.
2656/// \param[in] frame RooPlot to plot on
2657/// \param[in] fr The RooFitResult, where errors can be extracted
2658/// \param[in] Z The desired significance (width) of the error band
2659/// \param[in] params If non-zero, consider only the subset of the parameters in fr for the error evaluation
2660/// \param[in] argList Optional `RooCmdArg` that can be applied to a regular plotOn() operation
2661/// \param[in] linMethod By default (linMethod=true), a linearized error is shown.
2662/// \return The RooPlot the band was plotted on (for chaining of plotting commands).
2663///
2664/// The linearized error is calculated as follows:
2665/// \f[
2666/// \mathrm{error}(x) = Z * F_a(x) * \mathrm{Corr}(a,a') * F_{a'}^\mathrm{T}(x),
2667/// \f]
2668///
2669/// where
2670/// \f[
2671/// F_a(x) = \frac{ f(x,a+\mathrm{d}a) - f(x,a-\mathrm{d}a) }{2},
2672/// \f]
2673/// with \f$ f(x) \f$ the plotted curve and \f$ \mathrm{d}a \f$ taken from the fit result, and
2674/// \f$ \mathrm{Corr}(a,a') \f$ = the correlation matrix from the fit result, and \f$ Z \f$ = requested signifance (\f$ Z \sigma \f$ band)
2675///
2676/// The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters), but may
2677/// not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made
2678///
2679/// Alternatively, a more robust error is calculated using a sampling method. In this method a number of curves
2680/// is calculated with variations of the parameter values, as drawn from a multi-variate Gaussian p.d.f. that is constructed
2681/// from the fit results covariance matrix. The error(x) is determined by calculating a central interval that capture N% of the variations
2682/// 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
2683/// 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)
2684/// Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations, but may take (much)
2685/// longer to calculate.
2686
2687RooPlot* RooAbsReal::plotOnWithErrorBand(RooPlot* frame,const RooFitResult& fr, double Z,const RooArgSet* params, const RooLinkedList& argList, bool linMethod) const
2688{
2689 RooLinkedList plotArgListTmp(argList) ;
2690 RooCmdConfig::stripCmdList(plotArgListTmp,"VisualizeError,MoveToBack") ;
2691
2692 // Strip any 'internal normalization' arguments from list
2693 RooLinkedList plotArgList ;
2694 for (auto * cmd : static_range_cast<RooCmdArg*>(plotArgListTmp)) {
2695 if (std::string("Normalization")==cmd->GetName()) {
2696 if (((RooCmdArg*)cmd)->getInt(1)!=0) {
2697 } else {
2698 plotArgList.Add(cmd) ;
2699 }
2700 } else {
2701 plotArgList.Add(cmd) ;
2702 }
2703 }
2704
2705 // Function to plot a single curve, creating a copy of the plotArgList to
2706 // pass as plot command arguments. The "FillColor" command is removed because
2707 // it has no effect on plotting single curves and would cause a warning.
2708 auto plotFunc = [&](RooAbsReal const& absReal) {
2709 RooLinkedList tmp(plotArgList) ;
2710 RooCmdConfig::stripCmdList(tmp, "FillColor");
2711 absReal.plotOn(frame, tmp);
2712 };
2713
2714 // Generate central value curve
2715 plotFunc(*this);
2716 RooCurve* cenCurve = frame->getCurve() ;
2717 if(!cenCurve){
2718 coutE(Plotting) << ClassName() << "::" << GetName() << ":plotOnWithErrorBand: no curve for central value available" << std::endl;
2719 return frame;
2720 }
2721 frame->remove(nullptr,false) ;
2722
2723 RooCurve* band(nullptr) ;
2724 if (!linMethod) {
2725
2726 // *** Interval method ***
2727 //
2728 // Make N variations of parameters samples from V and visualize N% central interval where N% is defined from Z
2729
2730 // Clone self for internal use
2731 RooAbsReal* cloneFunc = (RooAbsReal*) cloneTree() ;
2732 RooArgSet cloneParams;
2733 cloneFunc->getObservables(&fr.floatParsFinal(), cloneParams) ;
2734 RooArgSet errorParams{cloneParams};
2735 if(params) {
2736 // clear and fill errorParams only with parameters that both in params and cloneParams
2737 cloneParams.selectCommon(*params, errorParams);
2738 }
2739
2740 // Generate 100 random parameter points distributed according to fit result covariance matrix
2741 RooAbsPdf* paramPdf = fr.createHessePdf(errorParams) ;
2742 Int_t n = Int_t(100./TMath::Erfc(Z/sqrt(2.))) ;
2743 if (n<100) n=100 ;
2744
2745 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") INFO: visualizing " << Z << "-sigma uncertainties in parameters "
2746 << errorParams << " from fit result " << fr.GetName() << " using " << n << " samplings." << std::endl ;
2747
2748 // Generate variation curves with above set of parameter values
2749 double ymin = frame->GetMinimum() ;
2750 double ymax = frame->GetMaximum() ;
2751 std::unique_ptr<RooDataSet> generatedData{paramPdf->generate(errorParams,n)};
2752 std::vector<RooCurve*> cvec ;
2753 for (int i=0 ; i<generatedData->numEntries() ; i++) {
2754 cloneParams.assign(*generatedData->get(i)) ;
2755 plotFunc(*cloneFunc);
2756 cvec.push_back(frame->getCurve()) ;
2757 frame->remove(nullptr,false) ;
2758 }
2759 frame->SetMinimum(ymin) ;
2760 frame->SetMaximum(ymax) ;
2761
2762
2763 // Generate upper and lower curve points from 68% interval around each point of central curve
2764 band = cenCurve->makeErrorBand(cvec,Z) ;
2765
2766 // Cleanup
2767 delete paramPdf ;
2768 delete cloneFunc ;
2769 for (std::vector<RooCurve*>::iterator i=cvec.begin() ; i!=cvec.end() ; ++i) {
2770 delete (*i) ;
2771 }
2772
2773 } else {
2774
2775 // *** Linear Method ***
2776 //
2777 // Make a one-sigma up- and down fluctation for each parameter and visualize
2778 // a from a linearized calculation as follows
2779 //
2780 // error(x) = F(a) C_aa' F(a')
2781 //
2782 // Where F(a) = (f(x,a+da) - f(x,a-da))/2
2783 // and C_aa' is the correlation matrix
2784
2785 // Strip out parameters with zero error
2786 RooArgList fpf_stripped;
2787 for (auto const* frv : static_range_cast<RooRealVar*>(fr.floatParsFinal())) {
2788 if (frv->getError() > frv->getVal() * std::numeric_limits<double>::epsilon()) {
2789 fpf_stripped.add(*frv);
2790 }
2791 }
2792
2793 // Clone self for internal use
2794 RooAbsReal* cloneFunc = (RooAbsReal*) cloneTree() ;
2795 RooArgSet cloneParams;
2796 cloneFunc->getObservables(&fpf_stripped, cloneParams) ;
2797 RooArgSet errorParams{cloneParams};
2798 if(params) {
2799 // clear and fill errorParams only with parameters that both in params and cloneParams
2800 cloneParams.selectCommon(*params, errorParams);
2801 }
2802
2803
2804 // Make list of parameter instances of cloneFunc in order of error matrix
2805 RooArgList paramList ;
2806 const RooArgList& fpf = fr.floatParsFinal() ;
2807 std::vector<int> fpf_idx ;
2808 for (Int_t i=0 ; i<fpf.getSize() ; i++) {
2809 RooAbsArg* par = errorParams.find(fpf[i].GetName()) ;
2810 if (par) {
2811 paramList.add(*par) ;
2812 fpf_idx.push_back(i) ;
2813 }
2814 }
2815
2816 std::vector<RooCurve*> plusVar, 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 (Int_t ivar=0 ; ivar<paramList.getSize() ; ivar++) {
2826
2827 RooRealVar& rrv = (RooRealVar&)fpf[fpf_idx[ivar]] ;
2828
2829 double cenVal = rrv.getVal() ;
2830 double errVal = sqrt(V(ivar,ivar)) ;
2831
2832 // Make Plus variation
2833 ((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 ((RooRealVar*)paramList.at(ivar))->setVal(cenVal-Z*errVal) ;
2843 plotFunc(*cloneFunc);
2844 minusVar.push_back(frame->getCurve()) ;
2845 frame->remove(nullptr,false) ;
2846
2847 ((RooRealVar*)paramList.at(ivar))->setVal(cenVal) ;
2848 }
2849
2850 TMatrixDSym C(paramList.getSize()) ;
2851 std::vector<double> errVec(paramList.size()) ;
2852 for (int i=0 ; i<paramList.getSize() ; i++) {
2853 errVec[i] = sqrt(V(i,i)) ;
2854 for (int j=i ; j<paramList.getSize() ; 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::Detail::owningPtr(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 = (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::Detail::owningPtr<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 && (matched.getSize()!=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 = (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::Detail::owningPtr(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 = (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 = (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 = (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::Detail::owningPtr(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 (int i=0 ; i<obs.getSize() ; 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 (int i=0 ; i<pars.getSize() ; 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.getSize()) {
4021 case 1: {
4022 RooRealVar* x = (RooRealVar*)obs.at(0) ;
4023 f = functor(obs,pars,nset) ;
4024 tf = new TF1(GetName(),f,x->getMin(),x->getMax(),pars.getSize()) ;
4025 break ;
4026 }
4027 case 2: {
4028 RooRealVar* x = (RooRealVar*)obs.at(0) ;
4029 RooRealVar* y = (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.getSize()) ;
4032 break ;
4033 }
4034 case 3: {
4035 RooRealVar* x = (RooRealVar*)obs.at(0) ;
4036 RooRealVar* y = (RooRealVar*)obs.at(1) ;
4037 RooRealVar* z = (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.getSize()) ;
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 (int i=0 ; i<pars.getSize() ; i++) {
4050 RooRealVar* p = (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////////////////////////////////////////////////////////////////////////////////
4197/// Calls RooAbsReal::createChi2(RooDataSet& data, const RooLinkedList& cmdList) and returns fit result.
4198///
4199/// For the list of possible commands in the `cmdList`, see
4200/// RooAbsReal::chi2FitTo(RooDataHist&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&) and
4201/// RooAbsPdf::chi2FitTo(RooDataHist&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&).
4202
4204{
4205 // Select the pdf-specific commands
4206 RooCmdConfig pc("RooAbsPdf::chi2FitTo(" + std::string(GetName()) + ")");
4207
4208 // Pull arguments to be passed to chi2 construction from list
4209 RooLinkedList fitCmdList(cmdList);
4210
4211 auto createChi2DataHistCmdArgs = "Range,RangeWithName,NumCPU,Optimize,IntegrateBins,ProjectedObservables,"
4212 "AddCoefRange,SplitRange,DataError,Extended";
4213 RooLinkedList chi2CmdList = pc.filterCmdList(fitCmdList, createChi2DataHistCmdArgs);
4214
4216
4217 // Process and check varargs
4218 pc.process(fitCmdList);
4219 if (!pc.ok(true)) {
4220 return nullptr;
4221 }
4222
4223 std::unique_ptr<RooAbsReal> chi2{createChi2(data, chi2CmdList)};
4225}
4226
4227
4228
4229
4230////////////////////////////////////////////////////////////////////////////////
4231/// Create a \f$ \chi^2 \f$ variable from a histogram and this function.
4232///
4233/// \param arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8 ordered arguments
4234///
4235/// The list of supported command arguments is given in the documentation for
4236/// 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&).
4237///
4238/// \param data Histogram with data
4239/// \return \f$ \chi^2 \f$ variable
4240
4242 const RooCmdArg &arg3, const RooCmdArg &arg4,
4243 const RooCmdArg &arg5, const RooCmdArg &arg6,
4244 const RooCmdArg &arg7, const RooCmdArg &arg8)
4245{
4246 // Construct Chi2
4248 std::string baseName = "chi2_" + std::string(GetName()) + "_" + data.GetName();
4249
4250 // Clear possible range attributes from previous fits.
4251 removeStringAttribute("fitrange");
4252
4253 auto chi2 = std::make_unique<RooChi2Var>(baseName.c_str(), baseName.c_str(), *this, data, arg1, arg2, arg3, arg4,
4254 arg5, arg6, arg7, arg8);
4256
4257 return RooFit::Detail::owningPtr(std::move(chi2));
4258}
4259
4260
4261
4262
4263////////////////////////////////////////////////////////////////////////////////
4264/// \see RooAbsReal::createChi2(RooDataHist&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4265/// \param data hist data
4266/// \param cmdList List with RooCmdArg() from the table
4267
4269{
4270 // Fill array of commands
4271 const RooCmdArg* cmds[8] ;
4272 Int_t i(0) ;
4273 for(auto * arg : static_range_cast<RooCmdArg*>(cmdList)) {
4274 cmds[i++] = arg ;
4275 }
4276 for (;i<8 ; i++) {
4277 cmds[i] = &RooCmdArg::none() ;
4278 }
4279
4280 return createChi2(data,*cmds[0],*cmds[1],*cmds[2],*cmds[3],*cmds[4],*cmds[5],*cmds[6],*cmds[7]) ;
4281
4282}
4283
4284
4285
4286
4287
4288////////////////////////////////////////////////////////////////////////////////
4289/// Perform a 2-D \f$ \chi^2 \f$ fit using a series of x and y values stored in the dataset `xydata`.
4290/// The y values can either be the event weights, or can be another column designated
4291/// by the YVar() argument. The y value must have errors defined for the \f$ \chi^2 \f$ to
4292/// be well defined.
4293///
4294/// <table>
4295/// <tr><th><th> Options to control construction of the chi-square
4296/// <tr><td> `YVar(RooRealVar& yvar)` <td> Designate given column in dataset as Y value
4297/// <tr><td> `Integrate(bool flag)` <td> Integrate function over range specified by X errors
4298/// rather than take value at bin center.
4299///
4300/// <tr><th><th> Options to control flow of fit procedure
4301/// <tr><td> `InitialHesse(bool flag)` <td> Flag controls if HESSE before MIGRAD as well, off by default
4302/// <tr><td> `Hesse(bool flag)` <td> Flag controls if HESSE is run after MIGRAD, on by default
4303/// <tr><td> `Minos(bool flag)` <td> Flag controls if MINOS is run after HESSE, on by default
4304/// <tr><td> `Minos(const RooArgSet& set)` <td> Only run MINOS on given subset of arguments
4305/// <tr><td> `Save(bool flag)` <td> Flag controls if RooFitResult object is produced and returned, off by default
4306/// <tr><td> `Strategy(Int_t flag)` <td> Set Minuit strategy (0 through 2, default is 1)
4307///
4308/// <tr><th><th> Options to control informational output
4309/// <tr><td> `Verbose(bool flag)` <td> Flag controls if verbose output is printed (NLL, parameter changes during fit
4310/// <tr><td> `Timer(bool flag)` <td> Time CPU and wall clock consumption of fit steps, off by default
4311/// <tr><td> `PrintLevel(Int_t level)` <td> Set Minuit print level (-1 through 3, default is 1). At -1 all RooFit informational
4312/// messages are suppressed as well
4313/// <tr><td> `Warnings(bool flag)` <td> Enable or disable MINUIT warnings (enabled by default)
4314/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per likelihood evaluation. A negative
4315/// value suppress output completely, a zero value will only print the error count per p.d.f component,
4316/// a positive value is will print details of each error up to numErr messages per p.d.f component.
4317/// </table>
4318
4320 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
4321 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
4322{
4324 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
4325 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
4326 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
4327 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
4328 return chi2FitTo(xydata,l) ;
4329}
4330
4331
4332
4333
4334////////////////////////////////////////////////////////////////////////////////
4335/// \copydoc RooAbsReal::chi2FitTo(RooDataSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4336
4338{
4339 // Select the pdf-specific commands
4340 RooCmdConfig pc("RooAbsPdf::chi2FitTo(" + std::string(GetName()) + ")");
4341
4342 // Pull arguments to be passed to chi2 construction from list
4343 RooLinkedList fitCmdList(cmdList);
4344 auto createChi2DataSetCmdArgs = "YVar,Integrate,RangeWithName,NumCPU,Verbose";
4345 RooLinkedList chi2CmdList = pc.filterCmdList(fitCmdList, createChi2DataSetCmdArgs);
4346
4348
4349 // Process and check varargs
4350 pc.process(fitCmdList);
4351 if (!pc.ok(true)) {
4352 return nullptr;
4353 }
4354
4355 std::unique_ptr<RooAbsReal> xychi2{createChi2(xydata, chi2CmdList)};
4356 return RooFit::Detail::owningPtr(RooFit::FitHelpers::minimize(*this, *xychi2, xydata, pc));
4357}
4358
4359
4360
4361
4362////////////////////////////////////////////////////////////////////////////////
4363/// Create a \f$ \chi^2 \f$ from a series of x and y values stored in a dataset.
4364/// The y values can either be the event weights (default), or can be another column designated
4365/// by the YVar() argument. The y value must have errors defined for the \f$ \chi^2 \f$ to
4366/// be well defined.
4367///
4368/// The following named arguments are supported
4369///
4370/// | | Options to control construction of the \f$ \chi^2 \f$
4371/// |-|-----------------------------------------
4372/// | `YVar(RooRealVar& yvar)` | Designate given column in dataset as Y value
4373/// | `Integrate(bool flag)` | Integrate function over range specified by X errors rather than take value at bin center.
4374///
4375
4377 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
4378 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
4379{
4381 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
4382 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
4383 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
4384 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
4385 return createChi2(data,l) ;
4386}
4387
4388
4389////////////////////////////////////////////////////////////////////////////////
4390/// See RooAbsReal::createChi2(RooDataSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
4391
4393{
4394 // Select the pdf-specific commands
4395 RooCmdConfig pc("RooAbsReal::createChi2(" + std::string(GetName()) + ")");
4396
4397 pc.defineInt("integrate", "Integrate", 0, 0);
4398 pc.defineObject("yvar", "YVar", 0, nullptr);
4399 pc.defineString("rangeName", "RangeWithName", 0, "", true);
4400 pc.defineInt("numcpu", "NumCPU", 0, 1);
4401 pc.defineInt("interleave", "NumCPU", 1, 0);
4402 pc.defineInt("verbose", "Verbose", 0, 0);
4403
4404 // Process and check varargs
4405 pc.process(cmdList);
4406 if (!pc.ok(true)) {
4407 return nullptr;
4408 }
4409
4410 // Decode command line arguments
4411 bool integrate = pc.getInt("integrate");
4412 RooRealVar *yvar = static_cast<RooRealVar *>(pc.getObject("yvar"));
4413 const char *rangeName = pc.getString("rangeName", 0, true);
4414 Int_t numcpu = pc.getInt("numcpu");
4415 Int_t numcpu_strategy = pc.getInt("interleave");
4416 // strategy 3 works only for RooSimultaneous.
4417 if (numcpu_strategy == 3 && !this->InheritsFrom("RooSimultaneous")) {
4418 coutW(Minimization) << "Cannot use a NumCpu Strategy = 3 when the pdf is not a RooSimultaneous, "
4419 "falling back to default strategy = 0"
4420 << endl;
4421 numcpu_strategy = 0;
4422 }
4423 RooFit::MPSplit interl = (RooFit::MPSplit)numcpu_strategy;
4424 bool verbose = pc.getInt("verbose");
4425
4427 cfg.rangeName = rangeName ? rangeName : "";
4428 cfg.nCPU = numcpu;
4429 cfg.interleave = interl;
4430 cfg.verbose = verbose;
4431 cfg.verbose = false;
4432
4433 std::string name = "chi2_" + std::string(GetName()) + "_" + data.GetName();
4434
4436 std::make_unique<RooXYChi2Var>(name.c_str(), name.c_str(), *this, data, yvar, integrate, cfg));
4437}
4438
4439
4440
4441////////////////////////////////////////////////////////////////////////////////
4442/// Return current evaluation error logging mode.
4443
4445{
4446 return _evalErrorMode ;
4447}
4448
4449////////////////////////////////////////////////////////////////////////////////
4450/// Set evaluation error logging mode. Options are
4451///
4452/// PrintErrors - Print each error through RooMsgService() as it occurs
4453/// CollectErrors - Accumulate errors, but do not print them. A subsequent call
4454/// to printEvalErrors() will print a summary
4455/// CountErrors - Accumulate error count, but do not print them.
4456///
4457
4459{
4460 _evalErrorMode = m;
4461}
4462
4463
4464////////////////////////////////////////////////////////////////////////////////
4465
4467{
4468 std::string plist ;
4469 for (auto const* arg : paramVars) {
4470 if (!dependsOnValue(*arg)) {
4471 coutW(InputArguments) << "RooAbsReal::setParameterizeIntegral(" << GetName()
4472 << ") function does not depend on listed parameter " << arg->GetName() << ", ignoring" << std::endl ;
4473 continue ;
4474 }
4475 if (!plist.empty()) plist += ":" ;
4476 plist += arg->GetName() ;
4477 }
4478 setStringAttribute("CACHEPARAMINT",plist.c_str()) ;
4479}
4480
4481
4482/** Base function for computing multiple values of a RooAbsReal
4483\param output The array where the results are stored
4484\param nEvents The number of events to be processed
4485\param dataMap A std::map containing the input data for the computations
4486**/
4487void RooAbsReal::computeBatch(double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const {
4488 // Find all servers that are serving real numbers to us, retrieve their batch data,
4489 // and switch them into "always clean" operating mode, so they return always the last-set value.
4490 struct ServerData {
4491 RooAbsArg* server;
4492 std::span<const double> batch;
4493 double oldValue;
4494 RooAbsArg::OperMode oldOperMode;
4495 bool oldValueDirty;
4496 bool oldShapeDirty;
4497 };
4498 std::vector<ServerData> ourServers;
4499 ourServers.reserve(servers().size());
4500
4501 for (auto server : servers()) {
4502 auto serverValues = dataMap.at(server);
4503 if(serverValues.empty()) continue;
4504
4505 // maybe we are still missing inhibit dirty here
4506 auto oldOperMode = server->operMode();
4507 // See note at the bottom of this function to learn why we can only set
4508 // the operation mode to "always clean" if there are no other value
4509 // clients.
4510 server->setOperMode(RooAbsArg::AClean);
4511 ourServers.push_back({server,
4512 serverValues,
4513 server->isCategory() ? static_cast<RooAbsCategory const*>(server)->getCurrentIndex() : static_cast<RooAbsReal const*>(server)->_value,
4514 oldOperMode,
4515 server->_valueDirty,
4516 server->_shapeDirty});
4517 // Prevent the server from evaluating; just return cached result, which we will side load:
4518 }
4519
4520
4521 // Make sure that we restore all state when we finish:
4522 struct RestoreStateRAII {
4523 RestoreStateRAII(std::vector<ServerData>& servers) :
4524 _servers{servers} { }
4525
4526 ~RestoreStateRAII() {
4527 for (auto& serverData : _servers) {
4528 serverData.server->setCachedValue(serverData.oldValue, true);
4529 serverData.server->setOperMode(serverData.oldOperMode);
4530 serverData.server->_valueDirty = serverData.oldValueDirty;
4531 serverData.server->_shapeDirty = serverData.oldShapeDirty;
4532 }
4533 }
4534
4535 std::vector<ServerData>& _servers;
4536 } restoreState{ourServers};
4537
4538
4539 // Advising to implement the batch interface makes only sense if the batch was not a scalar.
4540 // Otherwise, there would be no speedup benefit.
4542 coutI(FastEvaluations) << "The class " << ClassName() << " does not implement the faster batch evaluation interface."
4543 << " Consider requesting or implementing it to benefit from a speed up." << std::endl;
4544 }
4545
4546
4547 // For each event, write temporary values into our servers' caches, and run a single-value computation.
4548
4549 for (std::size_t i=0; i < nEvents; ++i) {
4550 for (auto& serv : ourServers) {
4551 serv.server->setCachedValue(serv.batch[std::min(i, serv.batch.size()-1)], false);
4552 }
4553
4554 output[i] = evaluate();
4555 }
4556}
4557
4558////////////////////////////////////////////////////////////////////////////////
4559/// This function defines the analytical integral translation for the class.
4560///
4561/// \param[in] code The code that decides the integrands.
4562/// \param[in] rangeName Name of the normalization range.
4563/// \param[in] ctx An object to manage auxiliary information for code-squashing.
4564///
4565/// \returns The representative code string of the integral for the given object.
4566std::string RooAbsReal::buildCallToAnalyticIntegral(Int_t /* code */, const char * /* rangeName */,
4567 RooFit::Detail::CodeSquashContext & /*ctx*/) const
4568{
4569 std::stringstream errorMsg;
4570 errorMsg << "An analytical integral function for class \"" << ClassName() << "\" has not yet been implemented.";
4571 coutE(Minimization) << errorMsg.str() << std::endl;
4572 throw std::runtime_error(errorMsg.str().c_str());
4573}
4574
4575double RooAbsReal::_DEBUG_getVal(const RooArgSet* normalisationSet) const {
4576
4577 const bool tmpFast = _fast;
4578 const double tmp = _value;
4579
4580 double fullEval = 0.;
4581 try {
4582 fullEval = getValV(normalisationSet);
4583 }
4584 catch (CachingError& error) {
4585 throw CachingError(std::move(error),
4586 FormatPdfTree() << *this);
4587 }
4588
4589 const double ret = (_fast && !_inhibitDirty) ? _value : fullEval;
4590
4591 if (std::isfinite(ret) && ( ret != 0. ? (ret - fullEval)/ret : ret - fullEval) > 1.E-9) {
4592#ifndef NDEBUG
4594#endif
4595 FormatPdfTree formatter;
4596 formatter << "--> (Scalar computation wrong here:)\n"
4597 << GetName() << " " << this << " _fast=" << tmpFast
4598 << "\n\tcached _value=" << std::setprecision(16) << tmp
4599 << "\n\treturning =" << ret
4600 << "\n\trecomputed =" << fullEval
4601 << "\n\tnew _value =" << _value << "] ";
4602 formatter << "\nServers:";
4603 for (const auto server : _serverList) {
4604 formatter << "\n ";
4605 server->printStream(formatter.stream(), kName | kClassName | kArgs | kExtras | kAddress | kValue, kInline);
4606 }
4607
4608 throw CachingError(formatter);
4609 }
4610
4611 return ret;
4612}
4613
4614
4615bool RooAbsReal::redirectServersHook(const RooAbsCollection & newServerList, bool mustReplaceAll,
4616 bool nameChange, bool isRecursiveStep)
4617{
4619 return RooAbsArg::redirectServersHook(newServerList, mustReplaceAll, nameChange, isRecursiveStep);
4620}
4621
4622
4623////////////////////////////////////////////////////////////////////////////////
4624
4626{
4627 for (RooAbsArg* arg : servers()) {
4628 if(auto realArg = dynamic_cast<RooAbsReal*>(arg)) {
4629 realArg->enableOffsetting(flag) ;
4630 }
4631 }
4632}
4633
4634
4635RooAbsReal::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:2468
R__EXTERN TSystem * gSystem
Definition TSystem.h:560
std::ostream & stream()
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:80
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:564
void setShapeDirty()
Notify that a shape-like property (e.g. binning) has changed.
Definition RooAbsArg.h:496
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:437
bool _fast
Definition RooAbsArg.h:718
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.
const RefCountList_t & servers() const
List of all servers of this object.
Definition RooAbsArg.h:209
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:109
void removeStringAttribute(const Text_t *key)
Delete a string attribute with a given key.
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition RooAbsArg.h:491
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:697
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:92
RefCountList_t _serverList
Definition RooAbsArg.h:635
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:252
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...
RooAbsBinning is the 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.
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
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.
RooAbsCollection is an 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.
Int_t getSize() const
Return the number of elements in the collection.
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
Definition