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