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