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