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