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