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