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