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