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