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