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