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