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