Logo ROOT  
Reference Guide
RooAbsPdf.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/** \class RooAbsPdf
19 \ingroup Roofitcore
20
21## RooAbsPdf, the base class of all PDFs
22
23RooAbsPdf is the abstract interface for all probability density
24functions. The class provides hybrid analytical/numerical
25normalization for its implementations, error tracing and a MC
26generator interface.
27
28### A Minimal PDF Implementation
29
30A minimal implementation of a PDF class derived from RooAbsPdf
31should override the `evaluate()` function. This function should
32return the PDF's value (which does not need to be normalised).
33
34
35#### Normalization/Integration
36
37Although the normalization of a PDF is an integral part of a
38probability density function, normalization is treated separately
39in RooAbsPdf. The reason is that a RooAbsPdf object is more than a
40PDF: it can be a building block for a more complex, composite PDF
41if any of its variables are functions instead of variables. In
42such cases the normalization of the composite may not be simply the
43integral over the dependents of the top level PDF as these are
44functions with potentially non-trivial Jacobian terms themselves.
45\note Therefore, no explicit attempt should be made to normalize the
46function output in evaluate(). In particular, normalisation constants
47can be omitted to speed up the function evaluations, and included later
48in the integration of the PDF (see below), which is called rarely in
49comparison to the `evaluate()` function.
50
51In addition, RooAbsPdf objects do not have a static concept of what
52variables are parameters and what variables are dependents (which
53need to be integrated over for a correct PDF normalization).
54Instead, the choice of normalization is always specified each time a
55normalized value is requested from the PDF via the getVal()
56method.
57
58RooAbsPdf manages the entire normalization logic of each PDF with
59help of a RooRealIntegral object, which coordinates the integration
60of a given choice of normalization. By default, RooRealIntegral will
61perform a fully numeric integration of all dependents. However,
62PDFs can advertise one or more (partial) analytical integrals of
63their function, and these will be used by RooRealIntegral, if it
64determines that this is safe (i.e. no hidden Jacobian terms,
65multiplication with other PDFs that have one or more dependents in
66commen etc).
67
68#### Implementing analytical integrals
69To implement analytical integrals, two functions must be implemented. First,
70
71```
72Int_t getAnalyticalIntegral(const RooArgSet& integSet, RooArgSet& anaIntSet)
73```
74should return the analytical integrals that are supported. `integSet`
75is the set of dependents for which integration is requested. The
76function should copy the subset of dependents it can analytically
77integrate to `anaIntSet`, and return a unique identification code for
78this integration configuration. If no integration can be
79performed, zero should be returned. Second,
80
81```
82Double_t analyticalIntegral(Int_t code)
83```
84
85implements the actual analytical integral(s) advertised by
86`getAnalyticalIntegral()`. This function will only be called with
87codes returned by `getAnalyticalIntegral()`, except code zero.
88
89The integration range for each dependent to be integrated can
90be obtained from the dependent's proxy functions `min()` and
91`max()`. Never call these proxy functions for any proxy not known to
92be a dependent via the integration code. Doing so may be
93ill-defined, e.g. in case the proxy holds a function, and will
94trigger an assert. Integrated category dependents should always be
95summed over all of their states.
96
97
98
99### Direct generation of observables
100
101Distributions for any PDF can be generated with the accept/reject method,
102but for certain PDFs, more efficient methods may be implemented. To
103implement direct generation of one or more observables, two
104functions need to be implemented, similar to those for analytical
105integrals:
106
107```
108Int_t getGenerator(const RooArgSet& generateVars, RooArgSet& directVars)
109```
110and
111```
112void generateEvent(Int_t code)
113```
114
115The first function advertises observables, for which distributions can be generated,
116similar to the way analytical integrals are advertised. The second
117function implements the actual generator for the advertised observables.
118
119The generated dependent values should be stored in the proxy
120objects. For this, the assignment operator can be used (i.e. `xProxy
121= 3.0` ). Never call assign to any proxy not known to be a dependent
122via the generation code. Doing so may be ill-defined, e.g. in case
123the proxy holds a function, and will trigger an assert.
124
125
126### Batched function evaluations (Advanced usage)
128To speed up computations with large numbers of data events in unbinned fits,
129it is beneficial to override `evaluateBatch()`. Like this, large batches of
130computations can be done, without having to call `evaluate()` for each single data event.
131`evaluateBatch()` should execute the same computation as `evaluate()`, but it
132may choose an implementation that is capable of SIMD computations.
133If evaluateBatch is not implemented, the classic and slower `evaluate()` will be
134called for each data event.
135*/
136
137#include "RooAbsPdf.h"
138
139#include "RooFit.h"
140#include "RooMsgService.h"
141#include "RooDataSet.h"
142#include "RooArgSet.h"
143#include "RooArgProxy.h"
144#include "RooRealProxy.h"
145#include "RooRealVar.h"
146#include "RooGenContext.h"
147#include "RooBinnedGenContext.h"
148#include "RooPlot.h"
149#include "RooCurve.h"
150#include "RooNLLVar.h"
151#include "RooMinuit.h"
152#include "RooCategory.h"
153#include "RooNameReg.h"
154#include "RooCmdConfig.h"
155#include "RooGlobalFunc.h"
156#include "RooAddition.h"
157#include "RooRandom.h"
158#include "RooNumIntConfig.h"
159#include "RooProjectedPdf.h"
160#include "RooInt.h"
161#include "RooCustomizer.h"
162#include "RooConstraintSum.h"
163#include "RooParamBinning.h"
164#include "RooNumCdf.h"
165#include "RooFitResult.h"
166#include "RooNumGenConfig.h"
167#include "RooCachedReal.h"
168#include "RooXYChi2Var.h"
169#include "RooChi2Var.h"
170#include "RooMinimizer.h"
171#include "RooRealIntegral.h"
172#include "RooWorkspace.h"
173
174#include "RooHelpers.h"
175#include "RooVDTHeaders.h"
176
177#include "TClass.h"
178#include "Riostream.h"
179#include "TMath.h"
180#include "TObjString.h"
181#include "TPaveText.h"
182#include "TList.h"
183#include "TH1.h"
184#include "TH2.h"
185#include "TMatrixD.h"
186#include "TMatrixDSym.h"
187#include "Math/CholeskyDecomp.h"
188#include "RooDerivative.h"
189
190#include <string>
191
192using namespace std;
193
195;
197;
198
201
202////////////////////////////////////////////////////////////////////////////////
203/// Default constructor
204
205RooAbsPdf::RooAbsPdf() : _norm(0), _normSet(0), _specGeneratorConfig(0)
206{
207 _errorCount = 0 ;
208 _negCount = 0 ;
209 _rawValue = 0 ;
211 _traceCount = 0 ;
212}
213
214
215
216////////////////////////////////////////////////////////////////////////////////
217/// Constructor with name and title only
218
219RooAbsPdf::RooAbsPdf(const char *name, const char *title) :
220 RooAbsReal(name,title), _norm(0), _normSet(0), _normMgr(this,10), _selectComp(kTRUE), _specGeneratorConfig(0)
221{
223 setTraceCounter(0) ;
224}
225
226
227
228////////////////////////////////////////////////////////////////////////////////
229/// Constructor with name, title, and plot range
230
231RooAbsPdf::RooAbsPdf(const char *name, const char *title,
232 Double_t plotMin, Double_t plotMax) :
233 RooAbsReal(name,title,plotMin,plotMax), _norm(0), _normSet(0), _normMgr(this,10), _selectComp(kTRUE), _specGeneratorConfig(0)
234{
236 setTraceCounter(0) ;
237}
238
239
240
241////////////////////////////////////////////////////////////////////////////////
242/// Copy constructor
243
244RooAbsPdf::RooAbsPdf(const RooAbsPdf& other, const char* name) :
245 RooAbsReal(other,name), _norm(0), _normSet(0),
246 _normMgr(other._normMgr,this), _selectComp(other._selectComp), _normRange(other._normRange)
247{
250
251 if (other._specGeneratorConfig) {
253 } else {
255 }
256}
257
258
259
260////////////////////////////////////////////////////////////////////////////////
261/// Destructor
262
264{
266}
267
268
269
270////////////////////////////////////////////////////////////////////////////////
271/// Return current value, normalized by integrating over
272/// the observables in `nset`. If `nset` is 0, the unnormalized value
273/// is returned. All elements of `nset` must be lvalues.
274///
275/// Unnormalized values are not cached.
276/// Doing so would be complicated as `_norm->getVal()` could
277/// spoil the cache and interfere with returning the cached
278/// return value. Since unnormalized calls are typically
279/// done in integration calls, there is no performance hit.
280
282{
283
284 // Special handling of case without normalization set (used in numeric integration of pdfs)
285 if (!nset) {
286 RooArgSet* tmp = _normSet ;
287 _normSet = 0 ;
288 Double_t val = evaluate() ;
289 _normSet = tmp ;
290 Bool_t error = traceEvalPdf(val) ;
291
292 if (error) {
293 return 0 ;
294 }
295 return val ;
296 }
297
298
299 // Process change in last data set used
300 Bool_t nsetChanged(kFALSE) ;
301 if (nset!=_normSet || _norm==0) {
302 nsetChanged = syncNormalization(nset) ;
303 }
304
305 // Return value of object. Calculated if dirty, otherwise cached value is returned.
306 if (isValueDirty() || nsetChanged || _norm->isValueDirty()) {
307
308 // Evaluate numerator
309 Double_t rawVal = evaluate() ;
310 Bool_t error = traceEvalPdf(rawVal) ; // Error checking and printing
311
312 // Evaluate denominator
313 Double_t normVal(_norm->getVal()) ;
314
315 if (normVal < 0. || (normVal == 0. && rawVal != 0)) {
316 //Unreasonable normalisations. A zero integral can be tolerated if the function vanishes.
317 error=kTRUE ;
318 std::stringstream msg;
319 msg << "p.d.f normalization integral is zero or negative: " << normVal;
320 logEvalError(msg.str().c_str());
321 }
322
323 // Raise global error flag if problems occur
324 if (error || (rawVal == 0. && normVal == 0.)) {
325 _value = 0 ;
326 } else {
327 _value = rawVal / normVal ;
328 }
329
331 }
332
333 return _value ;
334}
335
336
337////////////////////////////////////////////////////////////////////////////////
338/// Compute batch of values for given range, and normalise by integrating over
339/// the observables in `nset`. If `nset` is `nullptr`, unnormalized values
340/// are returned. All elements of `nset` must be lvalues.
341///
342/// \param[in] begin Begin of the batch.
343/// \param[in] maxSize Size of the requested range. May come out smaller.
344/// \param[in] normSet If not nullptr, normalise results by integrating over
345/// the variables in this set. The normalisation is only computed once, and applied
346/// to the full batch.
347///
348RooSpan<const double> RooAbsPdf::getValBatch(std::size_t begin, std::size_t maxSize,
349 const RooArgSet* normSet) const
350{
351 // Some PDFs do preprocessing here:
352 getValV(normSet);
353
356 _allBatchesDirty = false;
357 }
358
359 // Special handling of case without normalization set (used in numeric integration of pdfs)
360 if (!normSet) {
361 RooArgSet* tmp = _normSet ;
362 _normSet = nullptr;
363 auto outputs = evaluateBatch(begin, maxSize);
364 maxSize = outputs.size();
365 _normSet = tmp;
366
368
369 return outputs;
370 }
371
372
373 // Process change in last data set used
374 Bool_t nsetChanged(kFALSE) ;
375 if (normSet != _normSet || _norm == nullptr) {
376 nsetChanged = syncNormalization(normSet) ;
377 }
378
379 // TODO wait if batch is computing?
381 || nsetChanged || _norm->isValueDirty()) {
382
383 auto outputs = evaluateBatch(begin, maxSize);
384 maxSize = outputs.size();
385
386 // Evaluate denominator
387 const double normVal = _norm->getVal();
388
389 if (normVal < 0.
390 || (normVal == 0. && std::any_of(outputs.begin(), outputs.end(), [](double val){return val != 0;}))) {
391 logEvalError(Form("p.d.f normalization integral is zero or negative."
392 "\n\tInt(%s) = %f", GetName(), normVal));
393 }
394
395 if (normVal != 1. && normVal > 0.) {
396 const double invNorm = 1./normVal;
397 for (double& val : outputs) { //CHECK_VECTORISE
398 val *= invNorm;
399 }
400 }
401
403 }
404
405 const auto ret = _batchData.getBatch(begin, maxSize);
406
407 return ret;
408}
409
410////////////////////////////////////////////////////////////////////////////////
411/// Analytical integral with normalization (see RooAbsReal::analyticalIntegralWN() for further information)
412///
413/// This function applies the normalization specified by 'normSet' to the integral returned
414/// by RooAbsReal::analyticalIntegral(). The passthrough scenario (code=0) is also changed
415/// to return a normalized answer
416
417Double_t RooAbsPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
418{
419 cxcoutD(Eval) << "RooAbsPdf::analyticalIntegralWN(" << GetName() << ") code = " << code << " normset = " << (normSet?*normSet:RooArgSet()) << endl ;
420
421
422 if (code==0) return getVal(normSet) ;
423 if (normSet) {
424 return analyticalIntegral(code,rangeName) / getNorm(normSet) ;
425 } else {
426 return analyticalIntegral(code,rangeName) ;
427 }
428}
429
430
431
432////////////////////////////////////////////////////////////////////////////////
433/// Check that passed value is positive and not 'not-a-number'. If
434/// not, print an error, until the error counter reaches its set
435/// maximum.
436
438{
439 // check for a math error or negative value
440 Bool_t error(kFALSE) ;
441 if (TMath::IsNaN(value)) {
442 logEvalError(Form("p.d.f value is Not-a-Number (%f), forcing value to zero",value)) ;
443 error=kTRUE ;
444 }
445 if (value<0) {
446 logEvalError(Form("p.d.f value is less than zero (%f), forcing value to zero",value)) ;
447 error=kTRUE ;
448 }
449
450 // do nothing if we are no longer tracing evaluations and there was no error
451 if(!error) return error ;
452
453 // otherwise, print out this evaluations input values and result
454 if(++_errorCount <= 10) {
455 cxcoutD(Tracing) << "*** Evaluation Error " << _errorCount << " ";
456 if(_errorCount == 10) cxcoutD(Tracing) << "(no more will be printed) ";
457 }
458 else {
459 return error ;
460 }
461
462 Print() ;
463 return error ;
464}
465
466
467
468////////////////////////////////////////////////////////////////////////////////
469/// Return the integral of this PDF over all observables listed in 'nset'.
470
472{
473 if (!nset) return 1 ;
474
476 if (_verboseEval>1) cxcoutD(Tracing) << IsA()->GetName() << "::getNorm(" << GetName() << "): norm(" << _norm << ") = " << _norm->getVal() << endl ;
477
478 Double_t ret = _norm->getVal() ;
479 if (ret==0.) {
480 if(++_errorCount <= 10) {
481 coutW(Eval) << "RooAbsPdf::getNorm(" << GetName() << ":: WARNING normalization is zero, nset = " ; nset->Print("1") ;
482 if(_errorCount == 10) coutW(Eval) << "RooAbsPdf::getNorm(" << GetName() << ") INFO: no more messages will be printed " << endl ;
483 }
484 }
485
486 return ret ;
487}
488
489
490
491////////////////////////////////////////////////////////////////////////////////
492/// Return pointer to RooAbsReal object that implements calculation of integral over observables iset in range
493/// rangeName, optionally taking the integrand normalized over observables nset
494
495const RooAbsReal* RooAbsPdf::getNormObj(const RooArgSet* nset, const RooArgSet* iset, const TNamed* rangeName) const
496{
497
498 // Check normalization is already stored
499 CacheElem* cache = (CacheElem*) _normMgr.getObj(nset,iset,0,rangeName) ;
500 if (cache) {
501 return cache->_norm ;
502 }
503
504 // If not create it now
505 RooArgSet* depList = getObservables(iset) ;
506 RooAbsReal* norm = createIntegral(*depList,*nset, *getIntegratorConfig(), RooNameReg::str(rangeName)) ;
507 delete depList ;
508
509 // Store it in the cache
510 cache = new CacheElem(*norm) ;
511 _normMgr.setObj(nset,iset,cache,rangeName) ;
512
513 // And return the newly created integral
514 return norm ;
515}
516
517
518
519////////////////////////////////////////////////////////////////////////////////
520/// Verify that the normalization integral cached with this PDF
521/// is valid for given set of normalization observables.
522///
523/// If not, the cached normalization integral (if any) is deleted
524/// and a new integral is constructed for use with 'nset'.
525/// Elements in 'nset' can be discrete and real, but must be lvalues.
526///
527/// For functions that declare to be self-normalized by overloading the
528/// selfNormalized() function, a unit normalization is always constructed.
529
530Bool_t RooAbsPdf::syncNormalization(const RooArgSet* nset, Bool_t adjustProxies) const
531{
532
533// cout << IsA()->GetName() << "::syncNormalization(" << GetName() << ") nset = " << nset << " = " << (nset?*nset:RooArgSet()) << endl ;
534
535 _normSet = (RooArgSet*) nset ;
536
537 // Check if data sets are identical
538 CacheElem* cache = (CacheElem*) _normMgr.getObj(nset) ;
539 if (cache) {
540
541 Bool_t nsetChanged = (_norm!=cache->_norm) ;
542 _norm = cache->_norm ;
543
544
545// cout << "returning existing object " << _norm->GetName() << endl ;
546
547 if (nsetChanged && adjustProxies) {
548 // Update dataset pointers of proxies
549 ((RooAbsPdf*) this)->setProxyNormSet(nset) ;
550 }
551
552 return nsetChanged ;
553 }
554
555 // Update dataset pointers of proxies
556 if (adjustProxies) {
557 ((RooAbsPdf*) this)->setProxyNormSet(nset) ;
558 }
559
560 std::unique_ptr<RooArgSet> depList(getObservables(nset));
561
562 if (_verboseEval>0) {
563 if (!selfNormalized()) {
564 cxcoutD(Tracing) << IsA()->GetName() << "::syncNormalization(" << GetName()
565 << ") recreating normalization integral " << endl ;
566 if (depList) depList->printStream(ccoutD(Tracing),kName|kValue|kArgs,kSingleLine) ; else ccoutD(Tracing) << "<none>" << endl ;
567 } else {
568 cxcoutD(Tracing) << IsA()->GetName() << "::syncNormalization(" << GetName() << ") selfNormalized, creating unit norm" << endl;
569 }
570 }
571
572 // Destroy old normalization & create new
573 if (selfNormalized() || !dependsOn(*depList)) {
574 TString ntitle(GetTitle()) ; ntitle.Append(" Unit Normalization") ;
575 TString nname(GetName()) ; nname.Append("_UnitNorm") ;
576 _norm = new RooRealVar(nname.Data(),ntitle.Data(),1) ;
577 } else {
578 const char* nr = (_normRangeOverride.Length()>0 ? _normRangeOverride.Data() : (_normRange.Length()>0 ? _normRange.Data() : 0)) ;
579
580// cout << "RooAbsPdf::syncNormalization(" << GetName() << ") rangeName for normalization is " << (nr?nr:"<null>") << endl ;
581 RooAbsReal* normInt = createIntegral(*depList,*getIntegratorConfig(),nr) ;
582 normInt->getVal() ;
583// cout << "resulting normInt = " << normInt->GetName() << endl ;
584
585 const char* cacheParamsStr = getStringAttribute("CACHEPARAMINT") ;
586 if (cacheParamsStr && strlen(cacheParamsStr)) {
587
588 RooArgSet* intParams = normInt->getVariables() ;
589
590 RooNameSet cacheParamNames ;
591 cacheParamNames.setNameList(cacheParamsStr) ;
592 RooArgSet* cacheParams = cacheParamNames.select(*intParams) ;
593
594 if (cacheParams->getSize()>0) {
595 cxcoutD(Caching) << "RooAbsReal::createIntObj(" << GetName() << ") INFO: constructing " << cacheParams->getSize()
596 << "-dim value cache for integral over " << *depList << " as a function of " << *cacheParams << " in range " << (nr?nr:"<default>") << endl ;
597 string name = Form("%s_CACHE_[%s]",normInt->GetName(),cacheParams->contentsString().c_str()) ;
598 RooCachedReal* cachedIntegral = new RooCachedReal(name.c_str(),name.c_str(),*normInt,*cacheParams) ;
599 cachedIntegral->setInterpolationOrder(2) ;
600 cachedIntegral->addOwnedComponents(*normInt) ;
601 cachedIntegral->setCacheSource(kTRUE) ;
602 if (normInt->operMode()==ADirty) {
603 cachedIntegral->setOperMode(ADirty) ;
604 }
605 normInt= cachedIntegral ;
606 }
607
608 delete cacheParams ;
609 delete intParams ;
610 }
611 _norm = normInt ;
612 }
613
614 // Register new normalization with manager (takes ownership)
615 cache = new CacheElem(*_norm) ;
616 _normMgr.setObj(nset,cache) ;
617
618// cout << "making new object " << _norm->GetName() << endl ;
619
620 return kTRUE ;
621}
622
623
624
625////////////////////////////////////////////////////////////////////////////////
626/// WVE 08/21/01 Probably obsolete now.
627
628Bool_t RooAbsPdf::traceEvalHook(Double_t value) const
629{
630 // Floating point error checking and tracing for given float value
631
632 // check for a math error or negative value
633 Bool_t error= TMath::IsNaN(value) || (value < 0);
634
635 // do nothing if we are no longer tracing evaluations and there was no error
636 if(!error && _traceCount <= 0) return error ;
637
638 // otherwise, print out this evaluations input values and result
639 if(error && ++_errorCount <= 10) {
640 cxcoutD(Tracing) << "*** Evaluation Error " << _errorCount << " ";
641 if(_errorCount == 10) ccoutD(Tracing) << "(no more will be printed) ";
642 }
643 else if(_traceCount > 0) {
644 ccoutD(Tracing) << '[' << _traceCount-- << "] ";
645 }
646 else {
647 return error ;
648 }
649
650 Print() ;
651
652 return error ;
653}
654
655
656
657
658////////////////////////////////////////////////////////////////////////////////
659/// Reset error counter to given value, limiting the number
660/// of future error messages for this pdf to 'resetValue'
661
663{
664 _errorCount = resetValue ;
665 _negCount = resetValue ;
666}
667
668
669
670////////////////////////////////////////////////////////////////////////////////
671/// Reset trace counter to given value, limiting the
672/// number of future trace messages for this pdf to 'value'
673
675{
676 if (!allNodes) {
677 _traceCount = value ;
678 return ;
679 } else {
680 RooArgList branchList ;
681 branchNodeServerList(&branchList) ;
682 TIterator* iter = branchList.createIterator() ;
683 RooAbsArg* arg ;
684 while((arg=(RooAbsArg*)iter->Next())) {
685 RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
686 if (pdf) pdf->setTraceCounter(value,kFALSE) ;
687 }
688 delete iter ;
689 }
690
691}
692
693
694
695
696////////////////////////////////////////////////////////////////////////////////
697/// Return the log of the current value with given normalization
698/// An error message is printed if the argument of the log is negative.
699
701{
702 Double_t prob = getVal(nset) ;
703
704 if (fabs(prob)>1e6) {
705 coutW(Eval) << "RooAbsPdf::getLogVal(" << GetName() << ") WARNING: large likelihood value: " << prob << endl ;
706 }
707
708 if(prob < 0) {
709 logEvalError("getLogVal() top-level p.d.f evaluates to a negative number") ;
710
711 return std::numeric_limits<double>::quiet_NaN();
712 }
713
714 if(prob == 0) {
715 logEvalError("getLogVal() top-level p.d.f evaluates to zero") ;
716
717 return -std::numeric_limits<double>::infinity();
718 }
719
720 if (TMath::IsNaN(prob)) {
721 logEvalError("getLogVal() top-level p.d.f evaluates to NaN") ;
722
723 return -std::numeric_limits<double>::infinity();
724 }
725
726 return log(prob);
727}
728
729
730////////////////////////////////////////////////////////////////////////////////
731/// Check for infinity or NaN.
732/// \param[in] inputs Array to check
733/// \return True if either infinity or NaN were found.
734namespace {
735template<class T>
736bool checkInfNaNNeg(const T& inputs) {
737 // check for a math error or negative value
738 bool inf = false;
739 bool nan = false;
740 bool neg = false;
741
742 for (double val : inputs) { //CHECK_VECTORISE
743 inf |= !std::isfinite(val);
744 nan |= TMath::IsNaN(val); // Works also during fast math
745 neg |= val < 0;
746 }
747
748 return inf || nan || neg;
749}
750}
751
752
753////////////////////////////////////////////////////////////////////////////////
754/// Scan through outputs and fix+log all nans and negative values.
755/// \param[in/out] outputs Array to be scanned & fixed.
756/// \param[in] begin Begin of event range. Only needed to print the correct event number
757/// where the error occurred.
758void RooAbsPdf::logBatchComputationErrors(RooSpan<const double>& outputs, std::size_t begin) const {
759 for (unsigned int i=0; i<outputs.size(); ++i) {
760 const double value = outputs[i];
761 if (TMath::IsNaN(outputs[i])) {
762 logEvalError(Form("p.d.f value of (%s) is Not-a-Number (%f) for entry %zu",
763 GetName(), value, begin+i));
764 } else if (!std::isfinite(outputs[i])){
765 logEvalError(Form("p.d.f value of (%s) is (%f) for entry %zu",
766 GetName(), value, begin+i));
767 } else if (outputs[i] < 0.) {
768 logEvalError(Form("p.d.f value of (%s) is less than zero (%f) for entry %zu",
769 GetName(), value, begin+i));
770 }
771 }
772}
773
774////////////////////////////////////////////////////////////////////////////////
775/// Compute the log-likelihoods for all events in the requested batch.
776/// The arguments are passed over to getValBatch().
777/// \param[in] begin Start of the batch.
778/// \param[in] size Maximum size of the batch. Depending on data layout and memory, the batch
779/// may come back smaller.
780/// \return Returns a batch of doubles that contains the log probabilities.
781RooSpan<const double> RooAbsPdf::getLogValBatch(std::size_t begin, std::size_t maxSize,
782 const RooArgSet* normSet) const
783{
784 auto pdfValues = getValBatch(begin, maxSize, normSet);
785
786 if (checkInfNaNNeg(pdfValues)) {
787 logBatchComputationErrors(pdfValues, begin);
788 }
789
790 auto output = _batchData.makeWritableBatchUnInit(begin, pdfValues.size());
791
792 for (std::size_t i = 0; i < pdfValues.size(); ++i) { //CHECK_VECTORISE
793 const double prob = pdfValues[i];
794
795 double theLog = _rf_fast_log(prob);
796
797 if (prob < 0) {
798 theLog = std::numeric_limits<double>::quiet_NaN();
799 } else if (prob == 0 || TMath::IsNaN(prob)) {
800 theLog = -std::numeric_limits<double>::infinity();
801 }
802
803 output[i] = theLog;
804 }
805
806 return output;
807}
808
809////////////////////////////////////////////////////////////////////////////////
810/// Return the extended likelihood term (\f$ N_\mathrm{expect} - N_\mathrm{observed} \cdot \log(N_\mathrm{expect} \f$)
811/// of this PDF for the given number of observed events.
812///
813/// For successful operation, the PDF implementation must indicate that
814/// it is extendable by overloading `canBeExtended()`, and must
815/// implement the `expectedEvents()` function.
816
818{
819 // check if this PDF supports extended maximum likelihood fits
820 if(!canBeExtended()) {
821 coutE(InputArguments) << fName << ": this PDF does not support extended maximum likelihood"
822 << endl;
823 return 0;
824 }
825
826 Double_t expected= expectedEvents(nset);
827 if(expected < 0) {
828 coutE(InputArguments) << fName << ": calculated negative expected events: " << expected
829 << endl;
830 return 0;
831 }
832
833
834 // Explicitly handle case Nobs=Nexp=0
835 if (fabs(expected)<1e-10 && fabs(observed)<1e-10) {
836 return 0 ;
837 }
838
839 // Check for errors in Nexpected
840 if (expected<0 || TMath::IsNaN(expected)) {
841 logEvalError("extendedTerm #expected events is <0 or NaN") ;
842 return 0 ;
843 }
844
845 // calculate and return the negative log-likelihood of the Poisson
846 // factor for this dataset, dropping the constant log(observed!)
847 // Double_t extra=0;
848 // if(observed<1000000) {
849 // Double_t Delta1 = (expected-observed);
850 // Double_t Delta2 = observed*(log(expected)-log(observed+1));
851 // Double_t Const3 = 0.5*log(observed+1);
852 // extra= Delta1 - Delta2 + Const3;
853
854 // cout << " extra obs = " << observed << " exp = " << expected << endl ;
855 // cout << " extra orig = " << expected - observed*log(expected) << endl ;
856 // cout << " extra orig fix = " << expected - observed*log(expected) + log(TMath::Gamma(observed+1)) << endl ;
857 // cout << " extra new = " << extra << endl ;
858
859 // } else {
860 // Double_t sigma_square=expected;
861 // Double_t diff=observed-expected;
862 // extra=-log(sigma_square)/2 + (diff*diff)/(2*sigma_square);
863 // }
864
865 Double_t extra= expected - observed*log(expected);
866
867// cout << "RooAbsPdf::extendedTerm(" << GetName() << ") observed = " << observed << " expected = " << expected << endl ;
868
869 Bool_t trace(kFALSE) ;
870 if(trace) {
871 cxcoutD(Tracing) << fName << "::extendedTerm: expected " << expected << " events, got "
872 << observed << " events. extendedTerm = " << extra << endl;
873 }
874
875// cout << "RooAbsPdf::extendedTerm(" << GetName() << ") nExp = " << expected << " nObs = " << observed << endl ;
876 return extra;
877}
878
879
880
881////////////////////////////////////////////////////////////////////////////////
882/// Construct representation of -log(L) of PDF with given dataset. If dataset is unbinned, an unbinned likelihood is constructed. If the dataset
883/// is binned, a binned likelihood is constructed.
884///
885/// The following named arguments are supported
886///
887/// <table>
888/// <tr><th> Type of CmdArg <th> Effect on nll
889/// <tr><td> `ConditionalObservables(const RooArgSet& set)` <td> Do not normalize PDF over listed observables
890/// <tr><td> `Extended(Bool_t flag)` <td> Add extended likelihood term, off by default
891/// <tr><td> `Range(const char* name)` <td> Fit only data inside range with given name
892/// <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.
893/// Multiple comma separated range names can be specified.
894/// <tr><td> `SumCoefRange(const char* name)` <td> Set the range in which to interpret the coefficients of RooAddPdf components
895/// <tr><td> `NumCPU(int num, int strat)` <td> Parallelize NLL calculation on num CPUs
896/// <table>
897/// <tr><th> Strategy <th> Effect
898/// <tr><td> 0 = RooFit::BulkPartition (Default) <td> Divide events in N equal chunks
899/// <tr><td> 1 = RooFit::Interleave <td> Process event i%N in process N. Recommended for binned data with
900/// a substantial number of zero-bins, which will be distributed across processes more equitably in this strategy
901/// <tr><td> 2 = RooFit::SimComponents <td> Process each component likelihood of a RooSimultaneous fully in a single process
902/// and distribute components over processes. This approach can be benificial if normalization calculation time
903/// dominates the total computation time of a component (since the normalization calculation must be performed
904/// in each process in strategies 0 and 1. However beware that if the RooSimultaneous components do not share many
905/// parameters this strategy is inefficient: as most minuit-induced likelihood calculations involve changing
906/// a single parameter, only 1 of the N processes will be active most of the time if RooSimultaneous components
907/// do not share many parameters
908/// <tr><td> 3 = RooFit::Hybrid <td> Follow strategy 0 for all RooSimultaneous components, except those with less than
909/// 30 dataset entries, for which strategy 2 is followed.
910/// </table>
911/// <tr><td> `BatchMode(bool on)` <td> Batch evaluation mode. See createNLL().
912/// <tr><td> `Optimize(Bool_t flag)` <td> Activate constant term optimization (on by default)
913/// <tr><td> `SplitRange(Bool_t flag)` <td> Use separate fit ranges in a simultaneous fit. Actual range name for each subsample is assumed to
914/// be `rangeName_indexState`, where `indexState` is the state of the master index category of the simultaneous fit.
915/// Using `Range("range"), SplitRange()` as switches, different ranges could be set like this:
916/// ```
917/// myVariable.setRange("range_pi0", 135, 210);
918/// myVariable.setRange("range_gamma", 50, 210);
919/// ```
920/// <tr><td> `Constrain(const RooArgSet&pars)` <td> For p.d.f.s that contain internal parameter constraint terms (that is usually product PDFs, where one
921/// term of the product depends on parameters but not on the observable(s),), only apply constraints to the given subset of parameters.
922/// <tr><td> `ExternalConstraints(const RooArgSet& )` <td> Include given external constraints to likelihood by multiplying them with the original likelihood.
923/// <tr><td> `GlobalObservables(const RooArgSet&)` <td> Define the set of normalization observables to be used for the constraint terms.
924/// If none are specified the constrained parameters are used.
925/// <tr><td> `GlobalObservablesTag(const char* tagName)` <td> Define the set of normalization observables to be used for the constraint terms by
926/// a string attribute associated with pdf observables that match the given tagName.
927/// <tr><td> `Verbose(Bool_t flag)` <td> Controls RooFit informational messages in likelihood construction
928/// <tr><td> `CloneData(Bool flag)` <td> Use clone of dataset in NLL (default is true)
929/// <tr><td> `Offset(Bool_t)` <td> Offset likelihood by initial value (so that starting value of FCN in minuit is zero).
930/// This can improve numeric stability in simultaneous fits with components with large likelihood values
931/// </table>
932///
933///
934
935RooAbsReal* RooAbsPdf::createNLL(RooAbsData& data, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
936 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
937{
939 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
940 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
941 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
942 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
943 return createNLL(data,l) ;
944}
945
946
947
948
949////////////////////////////////////////////////////////////////////////////////
950/// Construct representation of -log(L) of PDFwith given dataset. If dataset is unbinned, an unbinned likelihood is constructed. If the dataset
951/// is binned, a binned likelihood is constructed.
952///
953/// See RooAbsPdf::createNLL(RooAbsData& data, RooCmdArg arg1, RooCmdArg arg2, RooCmdArg arg3, RooCmdArg arg4,
954/// RooCmdArg arg5, RooCmdArg arg6, RooCmdArg arg7, RooCmdArg arg8)
955/// for documentation of options
956
958{
959
960 // Select the pdf-specific commands
961 RooCmdConfig pc(Form("RooAbsPdf::createNLL(%s)",GetName())) ;
962
963 pc.defineString("rangeName","RangeWithName",0,"",kTRUE) ;
964 pc.defineString("addCoefRange","SumCoefRange",0,"") ;
965 pc.defineString("globstag","GlobalObservablesTag",0,"") ;
966 pc.defineDouble("rangeLo","Range",0,-999.) ;
967 pc.defineDouble("rangeHi","Range",1,-999.) ;
968 pc.defineInt("splitRange","SplitRange",0,0) ;
969 pc.defineInt("ext","Extended",0,2) ;
970 pc.defineInt("numcpu","NumCPU",0,1) ;
971 pc.defineInt("interleave","NumCPU",1,0) ;
972 pc.defineInt("verbose","Verbose",0,0) ;
973 pc.defineInt("optConst","Optimize",0,0) ;
974 pc.defineInt("cloneData","CloneData", 0, 2);
975 pc.defineSet("projDepSet","ProjectedObservables",0,0) ;
976 pc.defineSet("cPars","Constrain",0,0) ;
977 pc.defineSet("glObs","GlobalObservables",0,0) ;
978// pc.defineInt("constrAll","Constrained",0,0) ;
979 pc.defineInt("doOffset","OffsetLikelihood",0,0) ;
980 pc.defineSet("extCons","ExternalConstraints",0,0) ;
981 pc.defineInt("BatchMode", "BatchMode", 0, 0);
982 pc.defineMutex("Range","RangeWithName") ;
983// pc.defineMutex("Constrain","Constrained") ;
984 pc.defineMutex("GlobalObservables","GlobalObservablesTag") ;
985
986 // Process and check varargs
987 pc.process(cmdList) ;
988 if (!pc.ok(kTRUE)) {
989 return 0 ;
990 }
991
992 // Decode command line arguments
993 const char* rangeName = pc.getString("rangeName",0,kTRUE) ;
994 const char* addCoefRangeName = pc.getString("addCoefRange",0,kTRUE) ;
995 const char* globsTag = pc.getString("globstag",0,kTRUE) ;
996 Int_t ext = pc.getInt("ext") ;
997 Int_t numcpu = pc.getInt("numcpu") ;
998 RooFit::MPSplit interl = (RooFit::MPSplit) pc.getInt("interleave") ;
999
1000 Int_t splitr = pc.getInt("splitRange") ;
1001 Bool_t verbose = pc.getInt("verbose") ;
1002 Int_t optConst = pc.getInt("optConst") ;
1003 Int_t cloneData = pc.getInt("cloneData") ;
1004 Int_t doOffset = pc.getInt("doOffset") ;
1005
1006 // If no explicit cloneData command is specified, cloneData is set to true if optimization is activated
1007 if (cloneData==2) {
1008 cloneData = optConst ;
1009 }
1010
1011
1012 RooArgSet* cPars = pc.getSet("cPars") ;
1013 RooArgSet* glObs = pc.getSet("glObs") ;
1014 if (pc.hasProcessed("GlobalObservablesTag")) {
1015 if (glObs) delete glObs ;
1016 RooArgSet* allVars = getVariables() ;
1017 glObs = (RooArgSet*) allVars->selectByAttrib(globsTag,kTRUE) ;
1018 coutI(Minimization) << "User-defined specification of global observables definition with tag named '" << globsTag << "'" << endl ;
1019 delete allVars ;
1020 } else if (!pc.hasProcessed("GlobalObservables")) {
1021
1022 // Neither GlobalObservables nor GlobalObservablesTag has been processed - try if a default tag is defined in the head node
1023 // Check if head not specifies default global observable tag
1024 const char* defGlobObsTag = getStringAttribute("DefaultGlobalObservablesTag") ;
1025 if (defGlobObsTag) {
1026 coutI(Minimization) << "p.d.f. provides built-in specification of global observables definition with tag named '" << defGlobObsTag << "'" << endl ;
1027 if (glObs) delete glObs ;
1028 RooArgSet* allVars = getVariables() ;
1029 glObs = (RooArgSet*) allVars->selectByAttrib(defGlobObsTag,kTRUE) ;
1030 }
1031 }
1032
1033
1034 Bool_t doStripDisconnected=kFALSE ;
1035
1036 // If no explicit list of parameters to be constrained is specified apply default algorithm
1037 // All terms of RooProdPdfs that do not contain observables and share a parameters with one or more
1038 // terms that do contain observables are added as constraints.
1039 if (!cPars) {
1040 cPars = getParameters(data,kFALSE) ;
1041 doStripDisconnected=kTRUE ;
1042 }
1043 const RooArgSet* extCons = pc.getSet("extCons") ;
1044
1045 // Process automatic extended option
1046 if (ext==2) {
1047 ext = ((extendMode()==CanBeExtended || extendMode()==MustBeExtended)) ? 1 : 0 ;
1048 if (ext) {
1049 coutI(Minimization) << "p.d.f. provides expected number of events, including extended term in likelihood." << endl ;
1050 }
1051 }
1052
1053 if (pc.hasProcessed("Range")) {
1054 Double_t rangeLo = pc.getDouble("rangeLo") ;
1055 Double_t rangeHi = pc.getDouble("rangeHi") ;
1056
1057 // Create range with name 'fit' with above limits on all observables
1058 RooArgSet* obs = getObservables(&data) ;
1059 for (auto arg : *obs) {
1060 RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
1061 if (rrv) rrv->setRange("fit",rangeLo,rangeHi) ;
1062 }
1063 // Set range name to be fitted to "fit"
1064 rangeName = "fit" ;
1065 }
1066
1067 RooArgSet projDeps ;
1068 RooArgSet* tmp = pc.getSet("projDepSet") ;
1069 if (tmp) {
1070 projDeps.add(*tmp) ;
1071 }
1072
1073 // Construct NLL
1075 RooAbsReal* nll ;
1076 string baseName = Form("nll_%s_%s",GetName(),data.GetName()) ;
1077 if (!rangeName || strchr(rangeName,',')==0) {
1078 // Simple case: default range, or single restricted range
1079 //cout<<"FK: Data test 1: "<<data.sumEntries()<<endl;
1080
1081 auto theNLL = new RooNLLVar(baseName.c_str(),"-log(likelihood)",
1082 *this,data,projDeps,ext,rangeName,addCoefRangeName,numcpu,interl,
1083 verbose,splitr,cloneData);
1084 theNLL->batchMode(pc.getInt("BatchMode"));
1085 nll = theNLL;
1086 } else {
1087 // Composite case: multiple ranges
1088 RooArgList nllList ;
1089 auto tokens = RooHelpers::tokenise(rangeName, ",");
1090 for (const auto& token : tokens) {
1091 auto nllComp = new RooNLLVar(Form("%s_%s",baseName.c_str(),token.c_str()),"-log(likelihood)",
1092 *this,data,projDeps,ext,token.c_str(),addCoefRangeName,numcpu,interl,
1093 verbose,splitr,cloneData);
1094 nllComp->batchMode(pc.getInt("BatchMode"));
1095 nllList.add(*nllComp) ;
1096 }
1097 nll = new RooAddition(baseName.c_str(),"-log(likelihood)",nllList,kTRUE) ;
1098 }
1100
1101 // Collect internal and external constraint specifications
1102 RooArgSet allConstraints ;
1103
1104 if (_myws && _myws->set(Form("CACHE_CONSTR_OF_PDF_%s_FOR_OBS_%s", GetName(), RooNameSet(*data.get()).content()))) {
1105
1106 // retrieve from cache
1107 const RooArgSet *constr =
1108 _myws->set(Form("CACHE_CONSTR_OF_PDF_%s_FOR_OBS_%s", GetName(), RooNameSet(*data.get()).content()));
1109 coutI(Minimization) << "createNLL picked up cached consraints from workspace with " << constr->getSize()
1110 << " entries" << endl;
1111 allConstraints.add(*constr);
1112
1113 } else {
1114
1115 if (cPars && cPars->getSize() > 0) {
1116 RooArgSet *constraints = getAllConstraints(*data.get(), *cPars, doStripDisconnected);
1117 allConstraints.add(*constraints);
1118 delete constraints;
1119 }
1120 if (extCons) {
1121 allConstraints.add(*extCons);
1122 }
1123
1124 // write to cache
1125 if (_myws) {
1126 // cout << "createNLL: creating cache for allconstraints=" << allConstraints << endl ;
1127 coutI(Minimization) << "createNLL: caching constraint set under name "
1128 << Form("CONSTR_OF_PDF_%s_FOR_OBS_%s", GetName(), RooNameSet(*data.get()).content())
1129 << " with " << allConstraints.getSize() << " entries" << endl;
1131 Form("CACHE_CONSTR_OF_PDF_%s_FOR_OBS_%s", GetName(), RooNameSet(*data.get()).content()), allConstraints);
1132 }
1133 }
1134
1135 // Include constraints, if any, in likelihood
1136 RooAbsReal* nllCons(0) ;
1137 if (allConstraints.getSize()>0 && cPars) {
1138
1139 coutI(Minimization) << " Including the following constraint terms in minimization: " << allConstraints << endl ;
1140 if (glObs) {
1141 coutI(Minimization) << "The following global observables have been defined: " << *glObs << endl ;
1142 }
1143 nllCons = new RooConstraintSum(Form("%s_constr",baseName.c_str()),"nllCons",allConstraints,glObs ? *glObs : *cPars) ;
1144 nllCons->setOperMode(ADirty) ;
1145 RooAbsReal* orignll = nll ;
1146
1147 nll = new RooAddition(Form("%s_with_constr",baseName.c_str()),"nllWithCons",RooArgSet(*nll,*nllCons)) ;
1148 nll->addOwnedComponents(RooArgSet(*orignll,*nllCons)) ;
1149 }
1150
1151
1152 if (optConst) {
1154 }
1155
1156 if (doStripDisconnected) {
1157 delete cPars ;
1158 }
1159
1160 if (doOffset) {
1161 nll->enableOffsetting(kTRUE) ;
1162 }
1163
1164 return nll ;
1165}
1166
1167
1168
1169
1170
1171
1172////////////////////////////////////////////////////////////////////////////////
1173/// Fit PDF to given dataset. If dataset is unbinned, an unbinned maximum likelihood is performed. If the dataset
1174/// is binned, a binned maximum likelihood is performed. By default the fit is executed through the MINUIT
1175/// commands MIGRAD, HESSE in succession.
1176/// \param[in] data Data to fit the PDF to
1177/// \param[in] arg1 One or more arguments to control the behaviour of the fit
1178/// \return RooFitResult * with the fit status if option Save() is used, 0 otherwise.
1179///
1180/// The following named arguments are supported
1181///
1182/// <table>
1183/// <tr><th> Type of CmdArg <th> Options to control construction of -log(L)
1184/// <tr><td> `ConditionalObservables(const RooArgSet& set)` <td> Do not normalize PDF over listed observables
1185/// <tr><td> `Extended(Bool_t flag)` <td> Add extended likelihood term, off by default
1186/// <tr><td> `Range(const char* name)` <td> Fit only data inside range with given name. Multiple comma-separated range names can be specified.
1187/// <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.
1188/// <tr><td> `SumCoefRange(const char* name)` <td> Set the range in which to interpret the coefficients of RooAddPdf components
1189/// <tr><td> `NumCPU(int num, int strat)` <td> Parallelize NLL calculation on `num` CPUs
1190/// <table>
1191/// <tr><th> Strategy <th> Effect
1192/// <tr><td> 0 = RooFit::BulkPartition (Default) <td> Divide events in N equal chunks
1193/// <tr><td> 1 = RooFit::Interleave <td> Process event i%N in process N. Recommended for binned data with
1194/// a substantial number of zero-bins, which will be distributed across processes more equitably in this strategy
1195/// <tr><td> 2 = RooFit::SimComponents <td> Process each component likelihood of a RooSimultaneous fully in a single process
1196/// and distribute components over processes. This approach can be benificial if normalization calculation time
1197/// dominates the total computation time of a component (since the normalization calculation must be performed
1198/// in each process in strategies 0 and 1. However beware that if the RooSimultaneous components do not share many
1199/// parameters this strategy is inefficient: as most minuit-induced likelihood calculations involve changing
1200/// a single parameter, only 1 of the N processes will be active most of the time if RooSimultaneous components
1201/// do not share many parameters
1202/// <tr><td> 3 = RooFit::Hybrid <td> Follow strategy 0 for all RooSimultaneous components, except those with less than
1203/// 30 dataset entries, for which strategy 2 is followed.
1204/// </table>
1205/// <tr><td> `SplitRange(Bool_t flag)` <td> Use separate fit ranges in a simultaneous fit. Actual range name for each subsample is assumed
1206/// to by `rangeName_indexState` where indexState is the state of the master index category of the simultaneous fit.
1207/// Using `Range("range"), SplitRange()` as switches, different ranges could be set like this:
1208/// ```
1209/// myVariable.setRange("range_pi0", 135, 210);
1210/// myVariable.setRange("range_gamma", 50, 210);
1211/// ```
1212/// <tr><td> `Constrain(const RooArgSet&pars)` <td> For p.d.f.s that contain internal parameter constraint terms (that is usually product PDFs, where one
1213/// term of the product depends on parameters but not on the observable(s),), only apply constraints to the given subset of parameters.
1214/// <tr><td> `ExternalConstraints(const RooArgSet& )` <td> Include given external constraints to likelihood by multiplying them with the original likelihood.
1215/// <tr><td> `GlobalObservables(const RooArgSet&)` <td> Define the set of normalization observables to be used for the constraint terms.
1216/// If none are specified the constrained parameters are used.
1217/// <tr><td> `Offset(Bool_t)` <td> Offset likelihood by initial value (so that starting value of FCN in minuit is zero).
1218/// This can improve numeric stability in simultaneously fits with components with large likelihood values
1219/// <tr><td> `BatchMode(bool on)` <td> **Experimental** batch evaluation mode. This computes a batch of likelihood values at a time,
1220/// uses faster math functions and possibly auto vectorisation (this depends on the compiler flags).
1221/// Depending on hardware capabilities, the compiler flags and whether a batch evaluation function was
1222/// implemented for the PDFs of the model, likelihood computations are 2x to 10x faster.
1223/// The relative difference of the single log-likelihoods w.r.t. the legacy mode is usually better than 1.E-12,
1224/// and fit parameters usually agree to better than 1.E-6.
1225///
1226/// <tr><th><th> Options to control flow of fit procedure
1227/// <tr><td> `Minimizer(type,algo)` <td> Choose minimization package and algorithm to use. Default is MINUIT/MIGRAD through the RooMinimizer interface,
1228/// but others can be specified (through RooMinimizer interface). Select OldMinuit to use MINUIT through the old RooMinuit interface
1229/// <table>
1230/// <tr><th> Type <th> Algorithm
1231/// <tr><td> OldMinuit <td> migrad, simplex, minimize (=migrad+simplex), migradimproved (=migrad+improve)
1232/// <tr><td> Minuit <td> migrad, simplex, minimize (=migrad+simplex), migradimproved (=migrad+improve)
1233/// <tr><td> Minuit2 <td> migrad, simplex, minimize, scan
1234/// <tr><td> GSLMultiMin <td> conjugatefr, conjugatepr, bfgs, bfgs2, steepestdescent
1235/// <tr><td> GSLSimAn <td> -
1236/// </table>
1237///
1238/// <tr><td> `InitialHesse(Bool_t flag)` <td> Flag controls if HESSE before MIGRAD as well, off by default
1239/// <tr><td> `Optimize(Bool_t flag)` <td> Activate constant term optimization of test statistic during minimization (on by default)
1240/// <tr><td> `Hesse(Bool_t flag)` <td> Flag controls if HESSE is run after MIGRAD, on by default
1241/// <tr><td> `Minos(Bool_t flag)` <td> Flag controls if MINOS is run after HESSE, off by default
1242/// <tr><td> `Minos(const RooArgSet& set)` <td> Only run MINOS on given subset of arguments
1243/// <tr><td> `Save(Bool_t flag)` <td> Flag controls if RooFitResult object is produced and returned, off by default
1244/// <tr><td> `Strategy(Int_t flag)` <td> Set Minuit strategy (0 to 2, default is 1)
1245/// <tr><td> `FitOptions(const char* optStr)` <td> \deprecated Steer fit with classic options string (for backward compatibility).
1246/// \attention Use of this option excludes use of any of the new style steering options.
1247///
1248/// <tr><td> `SumW2Error(Bool_t flag)` <td> Apply correction to errors and covariance matrix.
1249/// This uses two covariance matrices, one with the weights, the other with squared weights,
1250/// to obtain the correct errors for weighted likelihood fits. If this option is activated, the
1251/// corrected covariance matrix is calculated as \f$ V_\mathrm{corr} = V C^{-1} V \f$, where \f$ V \f$ is the original
1252/// covariance matrix and \f$ C \f$ is the inverse of the covariance matrix calculated using the
1253/// squared weights. This allows to switch between two interpretations of errors:
1254/// <table>
1255/// <tr><th> SumW2Error <th> Interpretation
1256/// <tr><td> true <td> The errors reflect the uncertainty of the Monte Carlo simulation.
1257/// Use this if you want to know how much accuracy you can get from the available Monte Carlo statistics.
1258///
1259/// **Example**: Simulation with 1000 events, the average weight is 0.1.
1260/// The errors are as big as if one fitted to 1000 events.
1261/// <tr><td> false <td> The errors reflect the errors of a dataset, which is as big as the sum of weights.
1262/// Use this if you want to know what statistical errors you would get if you had a dataset with as many
1263/// events as the (weighted) Monte Carlo simulation represents.
1264///
1265/// **Example** (Data as above):
1266/// The errors are as big as if one fitted to 100 events.
1267/// </table>
1268///
1269/// <tr><td> `PrefitDataFraction(double fraction)`
1270/// <td> Runs a prefit on a small dataset of size fraction*(actual data size). This can speed up fits
1271/// by finding good starting values for the parameters for the actual fit.
1272/// \warning Prefitting may give bad results when used in binned analysis.
1273///
1274/// <tr><th><th> Options to control informational output
1275/// <tr><td> `Verbose(Bool_t flag)` <td> Flag controls if verbose output is printed (NLL, parameter changes during fit).
1276/// <tr><td> `Timer(Bool_t flag)` <td> Time CPU and wall clock consumption of fit steps, off by default.
1277/// <tr><td> `PrintLevel(Int_t level)` <td> Set Minuit print level (-1 to 3, default is 1). At -1 all RooFit informational messages are suppressed as well.
1278/// See RooMinimizer::PrintLevel for the meaning of the levels.
1279/// <tr><td> `Warnings(Bool_t flag)` <td> Enable or disable MINUIT warnings (enabled by default)
1280/// <tr><td> `PrintEvalErrors(Int_t numErr)` <td> Control number of p.d.f evaluation errors printed per likelihood evaluation.
1281/// A negative value suppresses output completely, a zero value will only print the error count per p.d.f component,
1282/// a positive value will print details of each error up to `numErr` messages per p.d.f component.
1283/// </table>
1284///
1285
1286RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3, const RooCmdArg& arg4,
1287 const RooCmdArg& arg5, const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
1288{
1290 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
1291 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
1292 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
1293 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
1294 return fitTo(data,l) ;
1295}
1296
1297
1298
1299////////////////////////////////////////////////////////////////////////////////
1300/// Fit PDF to given dataset. If dataset is unbinned, an unbinned maximum likelihood is performed. If the dataset
1301/// is binned, a binned maximum likelihood is performed. By default the fit is executed through the MINUIT
1302/// commands MIGRAD, HESSE and MINOS in succession.
1303///
1304/// See RooAbsPdf::fitTo(RooAbsData&,RooCmdArg&,RooCmdArg&,RooCmdArg&,RooCmdArg&,RooCmdArg&,RooCmdArg&,RooCmdArg&,RooCmdArg&)
1305///
1306/// for documentation of options
1307
1309{
1310 // Select the pdf-specific commands
1311 RooCmdConfig pc(Form("RooAbsPdf::fitTo(%s)",GetName())) ;
1312
1313 RooLinkedList fitCmdList(cmdList) ;
1314 RooLinkedList nllCmdList = pc.filterCmdList(fitCmdList,"ProjectedObservables,Extended,Range,"
1315 "RangeWithName,SumCoefRange,NumCPU,SplitRange,Constrained,Constrain,ExternalConstraints,"
1316 "CloneData,GlobalObservables,GlobalObservablesTag,OffsetLikelihood,BatchMode");
1317
1318 pc.defineDouble("prefit", "Prefit",0,0);
1319 pc.defineString("fitOpt","FitOptions",0,"") ;
1320 pc.defineInt("optConst","Optimize",0,2) ;
1321 pc.defineInt("verbose","Verbose",0,0) ;
1322 pc.defineInt("doSave","Save",0,0) ;
1323 pc.defineInt("doTimer","Timer",0,0) ;
1324 pc.defineInt("plevel","PrintLevel",0,1) ;
1325 pc.defineInt("strat","Strategy",0,1) ;
1326 pc.defineInt("initHesse","InitialHesse",0,0) ;
1327 pc.defineInt("hesse","Hesse",0,1) ;
1328 pc.defineInt("minos","Minos",0,0) ;
1329 pc.defineInt("ext","Extended",0,2) ;
1330 pc.defineInt("numcpu","NumCPU",0,1) ;
1331 pc.defineInt("numee","PrintEvalErrors",0,10) ;
1332 pc.defineInt("doEEWall","EvalErrorWall",0,1) ;
1333 pc.defineInt("doWarn","Warnings",0,1) ;
1334 pc.defineInt("doSumW2","SumW2Error",0,-1) ;
1335 pc.defineInt("doAsymptoticError","AsymptoticError",0,-1) ;
1336 pc.defineInt("doOffset","OffsetLikelihood",0,0) ;
1337 pc.defineString("mintype","Minimizer",0,"Minuit") ;
1338 pc.defineString("minalg","Minimizer",1,"minuit") ;
1339 pc.defineObject("minosSet","Minos",0,0) ;
1340 pc.defineSet("cPars","Constrain",0,0) ;
1341 pc.defineSet("extCons","ExternalConstraints",0,0) ;
1342 pc.defineMutex("FitOptions","Verbose") ;
1343 pc.defineMutex("FitOptions","Save") ;
1344 pc.defineMutex("FitOptions","Timer") ;
1345 pc.defineMutex("FitOptions","Strategy") ;
1346 pc.defineMutex("FitOptions","InitialHesse") ;
1347 pc.defineMutex("FitOptions","Hesse") ;
1348 pc.defineMutex("FitOptions","Minos") ;
1349 pc.defineMutex("Range","RangeWithName") ;
1350 pc.defineMutex("InitialHesse","Minimizer") ;
1351
1352 // Process and check varargs
1353 pc.process(fitCmdList) ;
1354 if (!pc.ok(kTRUE)) {
1355 return 0 ;
1356 }
1357
1358 // Decode command line arguments
1359 Double_t prefit = pc.getDouble("prefit");
1360 const char* fitOpt = pc.getString("fitOpt",0,kTRUE) ;
1361 Int_t optConst = pc.getInt("optConst") ;
1362 Int_t verbose = pc.getInt("verbose") ;
1363 Int_t doSave = pc.getInt("doSave") ;
1364 Int_t doTimer = pc.getInt("doTimer") ;
1365 Int_t plevel = pc.getInt("plevel") ;
1366 Int_t strat = pc.getInt("strat") ;
1367 Int_t initHesse= pc.getInt("initHesse") ;
1368 Int_t hesse = pc.getInt("hesse") ;
1369 Int_t minos = pc.getInt("minos") ;
1370 Int_t numee = pc.getInt("numee") ;
1371 Int_t doEEWall = pc.getInt("doEEWall") ;
1372 Int_t doWarn = pc.getInt("doWarn") ;
1373 Int_t doSumW2 = pc.getInt("doSumW2") ;
1374 Int_t doAsymptotic = pc.getInt("doAsymptoticError");
1375 const RooArgSet* minosSet = static_cast<RooArgSet*>(pc.getObject("minosSet")) ;
1376#ifdef __ROOFIT_NOROOMINIMIZER
1377 const char* minType =0 ;
1378#else
1379 const char* minType = pc.getString("mintype","Minuit") ;
1380 const char* minAlg = pc.getString("minalg","minuit") ;
1381#endif
1382
1383 // Determine if the dataset has weights
1384 Bool_t weightedData = data.isNonPoissonWeighted() ;
1385
1386 // Warn user that a method to determine parameter uncertainties should be provided if weighted data is offered
1387 if (weightedData && doSumW2==-1 && doAsymptotic==-1) {
1388 coutW(InputArguments) << "RooAbsPdf::fitTo(" << GetName() << ") WARNING: a likelihood fit is requested of what appears to be weighted data.\n"
1389 << " While the estimated values of the parameters will always be calculated taking the weights into account,\n"
1390 << " there are multiple ways to estimate the errors of the parameters. You are advised to make an \n"
1391 << " explicit choice for the error calculation:\n"
1392 << " - Either provide SumW2Error(true), to calculate a sum-of-weights-corrected HESSE error matrix\n"
1393 << " (error will be proportional to the number of events in MC).\n"
1394 << " - Or provide SumW2Error(false), to return errors from original HESSE error matrix\n"
1395 << " (which will be proportional to the sum of the weights, i.e., a dataset with <sum of weights> events).\n"
1396 << " - Or provide AsymptoticError(true), to use the asymptotically correct expression\n"
1397 << " (for details see https://arxiv.org/abs/1911.01303)."
1398 << endl ;
1399 }
1400
1401
1402 // Warn user that sum-of-weights correction does not apply to MINOS errrors
1403 if (doSumW2==1 && minos) {
1404 coutW(InputArguments) << "RooAbsPdf::fitTo(" << GetName() << ") WARNING: sum-of-weights correction does not apply to MINOS errors" << endl ;
1405 }
1406 if (doAsymptotic==1 && minos) {
1407 coutW(InputArguments) << "RooAbsPdf::fitTo(" << GetName() << ") WARNING: asymptotic correction does not apply to MINOS errors" << endl ;
1408 }
1409
1410 if (prefit != 0) {
1411 size_t nEvents = static_cast<size_t>(prefit*data.numEntries());
1412 if (prefit > 0.5 || nEvents < 100) {
1413 oocoutW(this,InputArguments) << "PrefitDataFraction should be in suitable range."
1414 << "With the current PrefitDataFraction=" << prefit
1415 << ", the number of events would be " << nEvents<< " out of "
1416 << data.numEntries() << ". Skipping prefit..." << endl;
1417 }
1418 else {
1419 size_t step = data.numEntries()/nEvents;
1420 RooArgSet tinyVars(*data.get());
1421 RooRealVar weight("weight","weight",1);
1422
1423 if (data.isWeighted()) tinyVars.add(weight);
1424
1425 RooDataSet tiny("tiny", "tiny", tinyVars,
1426 data.isWeighted() ? RooFit::WeightVar(weight) : RooCmdArg());
1427
1428 for (int i=0; i<data.numEntries(); i+=step)
1429 {
1430 const RooArgSet *event = data.get(i);
1431 tiny.add(*event, data.weight());
1432 }
1433 RooLinkedList tinyCmdList(cmdList) ;
1434 pc.filterCmdList(tinyCmdList,"Prefit,Hesse,Minos,Verbose,Save,Timer");
1435 RooCmdArg hesse_option = RooFit::Hesse(false);
1436 RooCmdArg print_option = RooFit::PrintLevel(-1);
1437
1438 tinyCmdList.Add(&hesse_option);
1439 tinyCmdList.Add(&print_option);
1440
1441 fitTo(tiny,tinyCmdList);
1442 }
1443 }
1444
1445 RooAbsReal* nll = createNLL(data,nllCmdList) ;
1446 RooFitResult *ret = 0 ;
1447
1448 //avoid setting both SumW2 and Asymptotic for uncertainty correction
1449 if (doSumW2==1 && doAsymptotic==1) {
1450 coutE(InputArguments) << "RooAbsPdf::fitTo(" << GetName() << ") ERROR: Cannot compute both asymptotically correct and SumW2 errors." << endl ;
1451 return ret;
1452 }
1453
1454 // Instantiate MINUIT
1455
1456 if (string(minType)!="OldMinuit") {
1457
1458#ifndef __ROOFIT_NOROOMINIMIZER
1459 RooMinimizer m(*nll) ;
1460
1461 m.setMinimizerType(minType) ;
1462
1463 m.setEvalErrorWall(doEEWall) ;
1464 if (doWarn==0) {
1465 // m.setNoWarn() ; WVE FIX THIS
1466 }
1467
1468 m.setPrintEvalErrors(numee) ;
1469 if (plevel!=1) {
1470 m.setPrintLevel(plevel) ;
1471 }
1472
1473 if (optConst) {
1474 // Activate constant term optimization
1475 m.optimizeConst(optConst) ;
1476 }
1477
1478 if (fitOpt) {
1479
1480 // Play fit options as historically defined
1481 ret = m.fit(fitOpt) ;
1482
1483 } else {
1484
1485 if (verbose) {
1486 // Activate verbose options
1487 m.setVerbose(1) ;
1488 }
1489 if (doTimer) {
1490 // Activate timer options
1491 m.setProfile(1) ;
1492 }
1493
1494 if (strat!=1) {
1495 // Modify fit strategy
1496 m.setStrategy(strat) ;
1497 }
1498
1499 if (initHesse) {
1500 // Initialize errors with hesse
1501 m.hesse() ;
1502 }
1503
1504 // Minimize using chosen algorithm
1505 m.minimize(minType,minAlg) ;
1506
1507 if (hesse) {
1508 // Evaluate errors with Hesse
1509 m.hesse() ;
1510 }
1511
1512 //asymptotically correct approach
1513 if (doAsymptotic==1 && m.getNPar()>0) {
1514 //Calculated corrected errors for weighted likelihood fits
1515 std::unique_ptr<RooFitResult> rw(m.save());
1516 //Weighted inverse Hessian matrix
1517 const TMatrixDSym& matV = rw->covarianceMatrix();
1518 coutI(Fitting) << "RooAbsPdf::fitTo(" << GetName() << ") Calculating covariance matrix according to the asymptotically correct approach. If you find this method useful please consider citing https://arxiv.org/abs/1911.01303." << endl;
1519
1520 //Initialise matrix containing first derivatives
1521 TMatrixDSym num(rw->floatParsFinal().getSize());
1522 for (int k=0; k<rw->floatParsFinal().getSize(); k++)
1523 for (int l=0; l<rw->floatParsFinal().getSize(); l++)
1524 num(k,l) = 0.0;
1525 RooArgSet* obs = getObservables(data);
1526 //Create derivative objects
1527 std::vector<std::unique_ptr<RooDerivative> > derivatives;
1528 const RooArgList& floated = rw->floatParsFinal();
1529 std::unique_ptr<RooArgSet> floatingparams( (RooArgSet*)getParameters(data)->selectByAttrib("Constant", false) );
1530 for (const auto paramresult : floated) {
1531 auto paraminternal = static_cast<RooRealVar*>(floatingparams->find(*paramresult));
1532 assert(floatingparams->find(*paramresult)->IsA() == RooRealVar::Class());
1533 std::unique_ptr<RooDerivative> deriv( derivative(*paraminternal, *obs, 1) );
1534 derivatives.push_back(std::move(deriv));
1535 }
1536
1537 //Loop over data
1538 for (int j=0; j<data.numEntries(); j++) {
1539 //Sets obs to current data point, this is where the pdf will be evaluated
1540 *obs = *data.get(j);
1541 //Determine first derivatives
1542 std::vector<Double_t> diffs(floated.getSize(), 0.0);
1543 for (int k=0; k < floated.getSize(); k++) {
1544 const auto paramresult = static_cast<RooRealVar*>(floated.at(k));
1545 auto paraminternal = static_cast<RooRealVar*>(floatingparams->find(*paramresult));
1546 //first derivative to parameter k at best estimate point for this measurement
1547 Double_t diff = derivatives.at(k)->getVal();
1548 //need to reset to best fit point after differentiation
1549 *paraminternal = paramresult->getVal();
1550 diffs.at(k) = diff;
1551 }
1552 //Fill numerator matrix
1553 Double_t prob = getVal(obs);
1554 for (int k=0; k<floated.getSize(); k++) {
1555 for (int l=0; l<floated.getSize(); l++) {
1556 num(k,l) += data.weight()*data.weight()*diffs.at(k)*diffs.at(l)/(prob*prob);
1557 }
1558 }
1559 }
1560 num.Similarity(matV);
1561
1562 //Propagate corrected errors to parameters objects
1563 m.applyCovarianceMatrix(num);
1564 }
1565
1566 if (doSumW2==1 && m.getNPar()>0) {
1567 // Make list of RooNLLVar components of FCN
1568 RooArgSet* comps = nll->getComponents();
1569 vector<RooNLLVar*> nllComponents;
1570 nllComponents.reserve(comps->getSize());
1571 TIterator* citer = comps->createIterator();
1572 RooAbsArg* arg;
1573 while ((arg=(RooAbsArg*)citer->Next())) {
1574 RooNLLVar* nllComp = dynamic_cast<RooNLLVar*>(arg);
1575 if (!nllComp) continue;
1576 nllComponents.push_back(nllComp);
1577 }
1578 delete citer;
1579 delete comps;
1580
1581 // Calculated corrected errors for weighted likelihood fits
1582 RooFitResult* rw = m.save();
1583 for (vector<RooNLLVar*>::iterator it = nllComponents.begin(); nllComponents.end() != it; ++it) {
1584 (*it)->applyWeightSquared(kTRUE);
1585 }
1586 coutI(Fitting) << "RooAbsPdf::fitTo(" << GetName() << ") Calculating sum-of-weights-squared correction matrix for covariance matrix" << endl ;
1587 m.hesse();
1588 RooFitResult* rw2 = m.save();
1589 for (vector<RooNLLVar*>::iterator it = nllComponents.begin(); nllComponents.end() != it; ++it) {
1590 (*it)->applyWeightSquared(kFALSE);
1591 }
1592
1593 // Apply correction matrix
1594 const TMatrixDSym& matV = rw->covarianceMatrix();
1595 TMatrixDSym matC = rw2->covarianceMatrix();
1597 CholeskyDecompGenDim<Double_t> decomp(matC.GetNrows(), matC);
1598 if (!decomp) {
1599 coutE(Fitting) << "RooAbsPdf::fitTo(" << GetName()
1600 << ") ERROR: Cannot apply sum-of-weights correction to covariance matrix: correction matrix calculated with weight-squared is singular" <<endl ;
1601 } else {
1602 // replace C by its inverse
1603 decomp.Invert(matC);
1604 // the class lies about the matrix being symmetric, so fill in the
1605 // part above the diagonal
1606 for (int i = 0; i < matC.GetNrows(); ++i)
1607 for (int j = 0; j < i; ++j) matC(j, i) = matC(i, j);
1608 matC.Similarity(matV);
1609 // C now contiains V C^-1 V
1610 // Propagate corrected errors to parameters objects
1611 m.applyCovarianceMatrix(matC);
1612 }
1613
1614 delete rw;
1615 delete rw2;
1616 }
1617
1618 if (minos) {
1619 // Evaluate errs with Minos
1620 if (minosSet) {
1621 m.minos(*minosSet) ;
1622 } else {
1623 m.minos() ;
1624 }
1625 }
1626
1627 // Optionally return fit result
1628 if (doSave) {
1629 string name = Form("fitresult_%s_%s",GetName(),data.GetName()) ;
1630 string title = Form("Result of fit of p.d.f. %s to dataset %s",GetName(),data.GetName()) ;
1631 ret = m.save(name.c_str(),title.c_str()) ;
1632 }
1633
1634 }
1635 if (optConst) {
1636 m.optimizeConst(0) ;
1637 }
1638
1639#endif
1640
1641 } else {
1642
1643 RooMinuit m(*nll) ;
1644
1645 m.setEvalErrorWall(doEEWall) ;
1646 if (doWarn==0) {
1647 m.setNoWarn() ;
1648 }
1649
1650 m.setPrintEvalErrors(numee) ;
1651 if (plevel!=1) {
1652 m.setPrintLevel(plevel) ;
1653 }
1654
1655 if (optConst) {
1656 // Activate constant term optimization
1657 m.optimizeConst(optConst) ;
1658 }
1659
1660 if (fitOpt) {
1661
1662 // Play fit options as historically defined
1663 ret = m.fit(fitOpt) ;
1664
1665 } else {
1666
1667 if (verbose) {
1668 // Activate verbose options
1669 m.setVerbose(1) ;
1670 }
1671 if (doTimer) {
1672 // Activate timer options
1673 m.setProfile(1) ;
1674 }
1675
1676 if (strat!=1) {
1677 // Modify fit strategy
1678 m.setStrategy(strat) ;
1679 }
1680
1681 if (initHesse) {
1682 // Initialize errors with hesse
1683 m.hesse() ;
1684 }
1685
1686 // Minimize using migrad
1687 m.migrad() ;
1688
1689 if (hesse) {
1690 // Evaluate errors with Hesse
1691 m.hesse() ;
1692 }
1693
1694 //asymptotically correct approach
1695 if (doAsymptotic==1 && m.getNPar()>0) {
1696 //Calculated corrected errors for weighted likelihood fits
1697 std::unique_ptr<RooFitResult> rw(m.save());
1698 //Weighted inverse Hessian matrix
1699 const TMatrixDSym& matV = rw->covarianceMatrix();
1700 coutI(Fitting) << "RooAbsPdf::fitTo(" << GetName() << ") Calculating covariance matrix according to the asymptotically correct approach. If you find this method useful please consider citing https://arxiv.org/abs/1911.01303." << endl;
1701
1702 //Initialise matrix containing first derivatives
1703 TMatrixDSym num(rw->floatParsFinal().getSize());
1704 for (int k=0; k<rw->floatParsFinal().getSize(); k++)
1705 for (int l=0; l<rw->floatParsFinal().getSize(); l++)
1706 num(k,l) = 0.0;
1707 RooArgSet* obs = getObservables(data);
1708 //Create derivative objects
1709 std::vector<std::unique_ptr<RooDerivative> > derivatives;
1710 const RooArgList& floated = rw->floatParsFinal();
1711 std::unique_ptr<RooArgSet> floatingparams( (RooArgSet*)getParameters(data)->selectByAttrib("Constant", false) );
1712 for (int k=0; k<floated.getSize(); k++) {
1713 RooRealVar* paramresult = (RooRealVar*)floated.at(k);
1714 RooRealVar* paraminternal = (RooRealVar*)floatingparams->find(paramresult->getTitle());
1715 std::unique_ptr<RooDerivative> deriv( derivative(*paraminternal, *obs, 1) );
1716 derivatives.push_back(std::move(deriv));
1717 }
1718
1719 //Loop over data
1720 for (int j=0; j<data.numEntries(); j++) {
1721 //Sets obs to current data point, this is where the pdf will be evaluated
1722 *obs = *data.get(j);
1723 //Determine first derivatives
1724 std::vector<Double_t> diffs(floated.getSize(), 0.0);
1725 for (int k=0; k<floated.getSize(); k++) {
1726 RooRealVar* paramresult = (RooRealVar*)floated.at(k);
1727 RooRealVar* paraminternal = (RooRealVar*)floatingparams->find(paramresult->getTitle());
1728 //First derivative to parameter k at best estimate point for this measurement
1729 Double_t diff = derivatives.at(k)->getVal();
1730 //Need to reset to best fit point after differentiation
1731 *paraminternal = paramresult->getVal();
1732 diffs.at(k) = diff;
1733 }
1734 //Fill numerator matrix
1735 Double_t prob = getVal(obs);
1736 for (int k=0; k<floated.getSize(); k++) {
1737 for (int l=0; l<floated.getSize(); l++) {
1738 num(k,l) += data.weight()*data.weight()*diffs.at(k)*diffs.at(l)/(prob*prob);
1739 }
1740 }
1741 }
1742 num.Similarity(matV);
1743
1744 //Propagate corrected errors to parameters objects
1745 m.applyCovarianceMatrix(num);
1746 }
1747
1748 if (doSumW2==1 && m.getNPar()>0) {
1749
1750 // Make list of RooNLLVar components of FCN
1751 list<RooNLLVar*> nllComponents ;
1752 RooArgSet* comps = nll->getComponents() ;
1753 RooAbsArg* arg ;
1754 TIterator* citer = comps->createIterator() ;
1755 while((arg=(RooAbsArg*)citer->Next())) {
1756 RooNLLVar* nllComp = dynamic_cast<RooNLLVar*>(arg) ;
1757 if (nllComp) {
1758 nllComponents.push_back(nllComp) ;
1759 }
1760 }
1761 delete citer ;
1762 delete comps ;
1763
1764 // Calculated corrected errors for weighted likelihood fits
1765 RooFitResult* rw = m.save() ;
1766 for (list<RooNLLVar*>::iterator iter1=nllComponents.begin() ; iter1!=nllComponents.end() ; ++iter1) {
1767 (*iter1)->applyWeightSquared(kTRUE) ;
1768 }
1769 coutI(Fitting) << "RooAbsPdf::fitTo(" << GetName() << ") Calculating sum-of-weights-squared correction matrix for covariance matrix" << endl ;
1770 m.hesse() ;
1771 RooFitResult* rw2 = m.save() ;
1772 for (list<RooNLLVar*>::iterator iter2=nllComponents.begin() ; iter2!=nllComponents.end() ; ++iter2) {
1773 (*iter2)->applyWeightSquared(kFALSE) ;
1774 }
1775
1776 // Apply correction matrix
1777 const TMatrixDSym& matV = rw->covarianceMatrix();
1778 TMatrixDSym matC = rw2->covarianceMatrix();
1780 CholeskyDecompGenDim<Double_t> decomp(matC.GetNrows(), matC);
1781 if (!decomp) {
1782 coutE(Fitting) << "RooAbsPdf::fitTo(" << GetName()
1783 << ") ERROR: Cannot apply sum-of-weights correction to covariance matrix: correction matrix calculated with weight-squared is singular" <<endl ;
1784 } else {
1785 // replace C by its inverse
1786 decomp.Invert(matC);
1787 // the class lies about the matrix being symmetric, so fill in the
1788 // part above the diagonal
1789 for (int i = 0; i < matC.GetNrows(); ++i)
1790 for (int j = 0; j < i; ++j) matC(j, i) = matC(i, j);
1791 matC.Similarity(matV);
1792 // C now contiains V C^-1 V
1793 // Propagate corrected errors to parameters objects
1794 m.applyCovarianceMatrix(matC);
1795 }
1796
1797 delete rw ;
1798 delete rw2 ;
1799 }
1800
1801 if (minos) {
1802 // Evaluate errs with Minos
1803 if (minosSet) {
1804 m.minos(*minosSet) ;
1805 } else {
1806 m.minos() ;
1807 }
1808 }
1809
1810 // Optionally return fit result
1811 if (doSave) {
1812 string name = Form("fitresult_%s_%s",GetName(),data.GetName()) ;
1813 string title = Form("Result of fit of p.d.f. %s to dataset %s",GetName(),data.GetName()) ;
1814 ret = m.save(name.c_str(),title.c_str()) ;
1815 }
1816
1817 }
1818
1819 if (optConst) {
1820 m.optimizeConst(0) ;
1821 }
1822
1823 }
1824
1825
1826
1827 // Cleanup
1828 delete nll ;
1829 return ret ;
1830}
1831
1832
1833
1834////////////////////////////////////////////////////////////////////////////////
1835/// Calls RooAbsPdf::createChi2(RooDataSet& data, const RooLinkedList& cmdList) and returns fit result.
1836
1838{
1839 // Select the pdf-specific commands
1840 RooCmdConfig pc(Form("RooAbsPdf::chi2FitTo(%s)",GetName())) ;
1841
1842 // Pull arguments to be passed to chi2 construction from list
1843 RooLinkedList fitCmdList(cmdList) ;
1844 RooLinkedList chi2CmdList = pc.filterCmdList(fitCmdList,"Range,RangeWithName,NumCPU,Optimize,ProjectedObservables,AddCoefRange,SplitRange,DataError,Extended") ;
1845
1846 RooAbsReal* chi2 = createChi2(data,chi2CmdList) ;
1847 RooFitResult* ret = chi2FitDriver(*chi2,fitCmdList) ;
1848
1849 // Cleanup
1850 delete chi2 ;
1851 return ret ;
1852}
1853
1854
1855
1856
1857////////////////////////////////////////////////////////////////////////////////
1858/// Create a \f$ \chi^2 \f$ from a histogram and this function.
1859///
1860/// Options to control construction of the \f$ \chi^2 \f$
1861/// ------------------------------------------
1862/// <table>
1863/// <tr><th> Type of CmdArg <th> Effect on \f$ \chi^2 \f$
1864/// <tr><td> `Extended()` <td> Use expected number of events of an extended p.d.f as normalization
1865/// <tr><td> `DataError()` <td> Choose between:
1866/// - Expected error [RooAbsData::Expected]
1867/// - Observed error (e.g. Sum-of-weights) [RooAbsData::SumW2]
1868/// - Poisson interval [RooAbsData::Poisson]
1869/// - Default: Expected error for unweighted data, Sum-of-weights for weighted data [RooAbsData::Auto]
1870/// <tr><td> `NumCPU()` <td> Activate parallel processing feature
1871/// <tr><td> `Range()` <td> Fit only selected region
1872/// <tr><td> `SumCoefRange()` <td> Set the range in which to interpret the coefficients of RooAddPdf components
1873/// <tr><td> `SplitRange()` <td> Fit ranges used in different categories get named after the category.
1874/// Using `Range("range"), SplitRange()` as switches, different ranges could be set like this:
1875/// ```
1876/// myVariable.setRange("range_pi0", 135, 210);
1877/// myVariable.setRange("range_gamma", 50, 210);
1878/// ```
1879/// <tr><td> `ConditionalObservables()` <td> Define projected observables
1880/// </table>
1881
1883 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
1884 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
1885{
1886 RooLinkedList cmdList ;
1887 cmdList.Add((TObject*)&arg1) ; cmdList.Add((TObject*)&arg2) ;
1888 cmdList.Add((TObject*)&arg3) ; cmdList.Add((TObject*)&arg4) ;
1889 cmdList.Add((TObject*)&arg5) ; cmdList.Add((TObject*)&arg6) ;
1890 cmdList.Add((TObject*)&arg7) ; cmdList.Add((TObject*)&arg8) ;
1891
1892 RooCmdConfig pc(Form("RooAbsPdf::createChi2(%s)",GetName())) ;
1893 pc.defineString("rangeName","RangeWithName",0,"",kTRUE) ;
1894 pc.allowUndefined(kTRUE) ;
1895 pc.process(cmdList) ;
1896 if (!pc.ok(kTRUE)) {
1897 return 0 ;
1898 }
1899 const char* rangeName = pc.getString("rangeName",0,kTRUE) ;
1900
1901 // Construct Chi2
1903 RooAbsReal* chi2 ;
1904 string baseName = Form("chi2_%s_%s",GetName(),data.GetName()) ;
1905
1906 if (!rangeName || strchr(rangeName,',')==0) {
1907 // Simple case: default range, or single restricted range
1908
1909 chi2 = new RooChi2Var(baseName.c_str(),baseName.c_str(),*this,data,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
1910
1911 } else {
1912
1913 // Find which argument is RangeWithName
1914 const RooCmdArg* rarg(0) ;
1915 string rcmd = "RangeWithName" ;
1916 if (arg1.GetName()==rcmd) rarg = &arg1 ;
1917 if (arg2.GetName()==rcmd) rarg = &arg2 ;
1918 if (arg3.GetName()==rcmd) rarg = &arg3 ;
1919 if (arg4.GetName()==rcmd) rarg = &arg4 ;
1920 if (arg5.GetName()==rcmd) rarg = &arg5 ;
1921 if (arg6.GetName()==rcmd) rarg = &arg6 ;
1922 if (arg7.GetName()==rcmd) rarg = &arg7 ;
1923 if (arg8.GetName()==rcmd) rarg = &arg8 ;
1924
1925 // Composite case: multiple ranges
1926 RooArgList chi2List ;
1927 const size_t bufSize = strlen(rangeName)+1;
1928 char* buf = new char[bufSize] ;
1929 strlcpy(buf,rangeName,bufSize) ;
1930 char* token = strtok(buf,",") ;
1931 while(token) {
1932 RooCmdArg subRangeCmd = RooFit::Range(token) ;
1933 // Construct chi2 while substituting original RangeWithName argument with subrange argument created above
1934 RooAbsReal* chi2Comp = new RooChi2Var(Form("%s_%s",baseName.c_str(),token),"chi^2",*this,data,
1935 &arg1==rarg?subRangeCmd:arg1,&arg2==rarg?subRangeCmd:arg2,
1936 &arg3==rarg?subRangeCmd:arg3,&arg4==rarg?subRangeCmd:arg4,
1937 &arg5==rarg?subRangeCmd:arg5,&arg6==rarg?subRangeCmd:arg6,
1938 &arg7==rarg?subRangeCmd:arg7,&arg8==rarg?subRangeCmd:arg8) ;
1939 chi2List.add(*chi2Comp) ;
1940 token = strtok(0,",") ;
1941 }
1942 delete[] buf ;
1943 chi2 = new RooAddition(baseName.c_str(),"chi^2",chi2List,kTRUE) ;
1944 }
1946
1947
1948 return chi2 ;
1949}
1950
1951
1952
1953
1954////////////////////////////////////////////////////////////////////////////////
1955/// Argument-list version of RooAbsPdf::createChi2()
1956
1958{
1959 // Select the pdf-specific commands
1960 RooCmdConfig pc(Form("RooAbsPdf::fitTo(%s)",GetName())) ;
1961
1962 pc.defineInt("integrate","Integrate",0,0) ;
1963 pc.defineObject("yvar","YVar",0,0) ;
1964
1965 // Process and check varargs
1966 pc.process(cmdList) ;
1967 if (!pc.ok(kTRUE)) {
1968 return 0 ;
1969 }
1970
1971 // Decode command line arguments
1972 Bool_t integrate = pc.getInt("integrate") ;
1973 RooRealVar* yvar = (RooRealVar*) pc.getObject("yvar") ;
1974
1975 string name = Form("chi2_%s_%s",GetName(),data.GetName()) ;
1976
1977 if (yvar) {
1978 return new RooXYChi2Var(name.c_str(),name.c_str(),*this,data,*yvar,integrate) ;
1979 } else {
1980 return new RooXYChi2Var(name.c_str(),name.c_str(),*this,data,integrate) ;
1981 }
1982}
1983
1984
1985
1986
1987////////////////////////////////////////////////////////////////////////////////
1988/// Print value of p.d.f, also print normalization integral that was last used, if any
1989
1990void RooAbsPdf::printValue(ostream& os) const
1991{
1992 getVal() ;
1993
1994 if (_norm) {
1995 os << evaluate() << "/" << _norm->getVal() ;
1996 } else {
1997 os << evaluate() ;
1998 }
1999}
2000
2001
2002
2003////////////////////////////////////////////////////////////////////////////////
2004/// Print multi line detailed information of this RooAbsPdf
2005
2006void RooAbsPdf::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
2007{
2009 os << indent << "--- RooAbsPdf ---" << endl;
2010 os << indent << "Cached value = " << _value << endl ;
2011 if (_norm) {
2012 os << indent << " Normalization integral: " << endl ;
2013 TString moreIndent(indent) ; moreIndent.Append(" ") ;
2015 }
2016}
2017
2018
2019
2020////////////////////////////////////////////////////////////////////////////////
2021/// Return a binned generator context
2022
2024{
2025 return new RooBinnedGenContext(*this,vars,0,0,verbose) ;
2026}
2027
2028
2029////////////////////////////////////////////////////////////////////////////////
2030/// Interface function to create a generator context from a p.d.f. This default
2031/// implementation returns a 'standard' context that works for any p.d.f
2032
2034 const RooArgSet* auxProto, Bool_t verbose) const
2035{
2036 return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
2037}
2038
2039
2040////////////////////////////////////////////////////////////////////////////////
2041
2042RooAbsGenContext* RooAbsPdf::autoGenContext(const RooArgSet &vars, const RooDataSet* prototype, const RooArgSet* auxProto,
2043 Bool_t verbose, Bool_t autoBinned, const char* binnedTag) const
2044{
2045 if (prototype || (auxProto && auxProto->getSize()>0)) {
2046 return genContext(vars,prototype,auxProto,verbose);
2047 }
2048
2049 RooAbsGenContext *context(0) ;
2050 if ( (autoBinned && isBinnedDistribution(vars)) || ( binnedTag && strlen(binnedTag) && (getAttribute(binnedTag)||string(binnedTag)=="*"))) {
2051 context = binnedGenContext(vars,verbose) ;
2052 } else {
2053 context= genContext(vars,0,0,verbose);
2054 }
2055 return context ;
2056}
2057
2058
2059
2060////////////////////////////////////////////////////////////////////////////////
2061/// Generate a new dataset containing the specified variables with events sampled from our distribution.
2062/// Generate the specified number of events or expectedEvents() if not specified.
2063/// \param[in] whatVars Choose variables in which to generate events. Variables not listed here will remain
2064/// constant and not be used for event generation.
2065/// \param[in] argxx Optional RooCmdArg() to change behaviour of generate().
2066/// \return RooDataSet *, owned by caller.
2067///
2068/// Any variables of this PDF that are not in whatVars will use their
2069/// current values and be treated as fixed parameters. Returns zero
2070/// in case of an error.
2071///
2072/// <table>
2073/// <tr><th> Type of CmdArg <th> Effect on generate
2074/// <tr><td> `Name(const char* name)` <td> Name of the output dataset
2075/// <tr><td> `Verbose(Bool_t flag)` <td> Print informational messages during event generation
2076/// <tr><td> `NumEvent(int nevt)` <td> Generate specified number of events
2077/// <tr><td> `Extended()` <td> If no number of events to be generated is given,
2078/// use expected number of events from extended likelihood term.
2079/// This evidently only works for extended PDFs.
2080/// <tr><td> `GenBinned(const char* tag)` <td> Use binned generation for all component pdfs that have 'setAttribute(tag)' set
2081/// <tr><td> `AutoBinned(Bool_t flag)` <td> Automatically deploy binned generation for binned distributions (e.g. RooHistPdf, sums and products of
2082/// RooHistPdfs etc)
2083/// \note Datasets that are generated in binned mode are returned as weighted unbinned datasets. This means that
2084/// for each bin, there will be one event in the dataset with a weight corresponding to the (possibly randomised) bin content.
2085///
2086///
2087/// <tr><td> `AllBinned()` <td> As above, but for all components.
2088/// \note The notion of components is only meaningful for simultaneous PDFs
2089/// as binned generation is always executed at the top-level node for a regular
2090/// PDF, so for those it only mattes that the top-level node is tagged.
2091///
2092/// <tr><td> ProtoData(const RooDataSet& data, Bool_t randOrder)
2093/// <td> Use specified dataset as prototype dataset. If randOrder in ProtoData() is set to true,
2094/// the order of the events in the dataset will be read in a random order if the requested
2095/// number of events to be generated does not match the number of events in the prototype dataset.
2096/// \note If ProtoData() is used, the specified existing dataset as a prototype: the new dataset will contain
2097/// the same number of events as the prototype (unless otherwise specified), and any prototype variables not in
2098/// whatVars will be copied into the new dataset for each generated event and also used to set our PDF parameters.
2099/// The user can specify a number of events to generate that will override the default. The result is a
2100/// copy of the prototype dataset with only variables in whatVars randomized. Variables in whatVars that
2101/// are not in the prototype will be added as new columns to the generated dataset.
2102///
2103/// </table>
2104///
2105
2106RooDataSet *RooAbsPdf::generate(const RooArgSet& whatVars, const RooCmdArg& arg1,const RooCmdArg& arg2,
2107 const RooCmdArg& arg3,const RooCmdArg& arg4, const RooCmdArg& arg5,const RooCmdArg& arg6)
2108{
2109 // Select the pdf-specific commands
2110 RooCmdConfig pc(Form("RooAbsPdf::generate(%s)",GetName())) ;
2111 pc.defineObject("proto","PrototypeData",0,0) ;
2112 pc.defineString("dsetName","Name",0,"") ;
2113 pc.defineInt("randProto","PrototypeData",0,0) ;
2114 pc.defineInt("resampleProto","PrototypeData",1,0) ;
2115 pc.defineInt("verbose","Verbose",0,0) ;
2116 pc.defineInt("extended","Extended",0,0) ;
2117 pc.defineInt("nEvents","NumEvents",0,0) ;
2118 pc.defineInt("autoBinned","AutoBinned",0,1) ;
2119 pc.defineInt("expectedData","ExpectedData",0,0) ;
2120 pc.defineDouble("nEventsD","NumEventsD",0,-1.) ;
2121 pc.defineString("binnedTag","GenBinned",0,"") ;
2122 pc.defineMutex("GenBinned","ProtoData") ;
2123 pc.defineMutex("Extended", "NumEvents");
2124
2125 // Process and check varargs
2126 pc.process(arg1,arg2,arg3,arg4,arg5,arg6) ;
2127 if (!pc.ok(kTRUE)) {
2128 return 0 ;
2129 }
2130
2131 // Decode command line arguments
2132 RooDataSet* protoData = static_cast<RooDataSet*>(pc.getObject("proto",0)) ;
2133 const char* dsetName = pc.getString("dsetName") ;
2134 Bool_t verbose = pc.getInt("verbose") ;
2135 Bool_t randProto = pc.getInt("randProto") ;
2136 Bool_t resampleProto = pc.getInt("resampleProto") ;
2137 Bool_t extended = pc.getInt("extended") ;
2138 Bool_t autoBinned = pc.getInt("autoBinned") ;
2139 const char* binnedTag = pc.getString("binnedTag") ;
2140 Int_t nEventsI = pc.getInt("nEvents") ;
2141 Double_t nEventsD = pc.getInt("nEventsD") ;
2142 //Bool_t verbose = pc.getInt("verbose") ;
2143 Bool_t expectedData = pc.getInt("expectedData") ;
2144
2145 Double_t nEvents = (nEventsD>0) ? nEventsD : Double_t(nEventsI);
2146
2147 // Force binned mode for expected data mode
2148 if (expectedData) {
2149 binnedTag="*" ;
2150 }
2151
2152 if (extended) {
2153 if (nEvents == 0) nEvents = expectedEvents(&whatVars);
2154 } else if (nEvents==0) {
2155 cxcoutI(Generation) << "No number of events specified , number of events generated is "
2156 << GetName() << "::expectedEvents() = " << expectedEvents(&whatVars)<< endl ;
2157 }
2158
2159 if (extended && protoData && !randProto) {
2160 cxcoutI(Generation) << "WARNING Using generator option Extended() (Poisson distribution of #events) together "
2161 << "with a prototype dataset implies incomplete sampling or oversampling of proto data. "
2162 << "Set randomize flag in ProtoData() option to randomize prototype dataset order and thus "
2163 << "to randomize the set of over/undersampled prototype events for each generation cycle." << endl ;
2164 }
2165
2166
2167 // Forward to appropriate implementation
2168 RooDataSet* data ;
2169 if (protoData) {
2170 data = generate(whatVars,*protoData,Int_t(nEvents),verbose,randProto,resampleProto) ;
2171 } else {
2172 data = generate(whatVars,nEvents,verbose,autoBinned,binnedTag,expectedData, extended) ;
2173 }
2174
2175 // Rename dataset to given name if supplied
2176 if (dsetName && strlen(dsetName)>0) {
2177 data->SetName(dsetName) ;
2178 }
2179
2180 return data ;
2181}
2182
2183
2184
2185
2186
2187
2188////////////////////////////////////////////////////////////////////////////////
2189/// \note This method does not perform any generation. To generate according to generations specification call RooAbsPdf::generate(RooAbsPdf::GenSpec&) const
2190///
2191/// Details copied from RooAbsPdf::generate():
2192/// --------------------------------------------
2193/// \copydetails RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
2194
2196 const RooCmdArg& arg1,const RooCmdArg& arg2,
2197 const RooCmdArg& arg3,const RooCmdArg& arg4,
2198 const RooCmdArg& arg5,const RooCmdArg& arg6)
2199{
2200
2201 // Select the pdf-specific commands
2202 RooCmdConfig pc(Form("RooAbsPdf::generate(%s)",GetName())) ;
2203 pc.defineObject("proto","PrototypeData",0,0) ;
2204 pc.defineString("dsetName","Name",0,"") ;
2205 pc.defineInt("randProto","PrototypeData",0,0) ;
2206 pc.defineInt("resampleProto","PrototypeData",1,0) ;
2207 pc.defineInt("verbose","Verbose",0,0) ;
2208 pc.defineInt("extended","Extended",0,0) ;
2209 pc.defineInt("nEvents","NumEvents",0,0) ;
2210 pc.defineInt("autoBinned","AutoBinned",0,1) ;
2211 pc.defineString("binnedTag","GenBinned",0,"") ;
2212 pc.defineMutex("GenBinned","ProtoData") ;
2213
2214
2215 // Process and check varargs
2216 pc.process(arg1,arg2,arg3,arg4,arg5,arg6) ;
2217 if (!pc.ok(kTRUE)) {
2218 return 0 ;
2219 }
2220
2221 // Decode command line arguments
2222 RooDataSet* protoData = static_cast<RooDataSet*>(pc.getObject("proto",0)) ;
2223 const char* dsetName = pc.getString("dsetName") ;
2224 Int_t nEvents = pc.getInt("nEvents") ;
2225 Bool_t verbose = pc.getInt("verbose") ;
2226 Bool_t randProto = pc.getInt("randProto") ;
2227 Bool_t resampleProto = pc.getInt("resampleProto") ;
2228 Bool_t extended = pc.getInt("extended") ;
2229 Bool_t autoBinned = pc.getInt("autoBinned") ;
2230 const char* binnedTag = pc.getString("binnedTag") ;
2231
2232 RooAbsGenContext* cx = autoGenContext(whatVars,protoData,0,verbose,autoBinned,binnedTag) ;
2233
2234 return new GenSpec(cx,whatVars,protoData,nEvents,extended,randProto,resampleProto,dsetName) ;
2235}
2236
2237
2238////////////////////////////////////////////////////////////////////////////////
2239/// If many identical generation requests
2240/// are needed, e.g. in toy MC studies, it is more efficient to use the prepareMultiGen()/generate()
2241/// combination than calling the standard generate() multiple times as
2242/// initialization overhead is only incurred once.
2243
2245{
2246 //Int_t nEvt = spec._extended ? RooRandom::randomGenerator()->Poisson(spec._nGen) : spec._nGen ;
2247 //Int_t nEvt = spec._extended ? RooRandom::randomGenerator()->Poisson(spec._nGen==0?expectedEvents(spec._whatVars):spec._nGen) : spec._nGen ;
2248 //Int_t nEvt = spec._nGen == 0 ? RooRandom::randomGenerator()->Poisson(expectedEvents(spec._whatVars)) : spec._nGen;
2249
2250 Double_t nEvt = spec._nGen == 0 ? expectedEvents(spec._whatVars) : spec._nGen;
2251
2252 RooDataSet* ret = generate(*spec._genContext,spec._whatVars,spec._protoData, nEvt,kFALSE,spec._randProto,spec._resampleProto,
2253 spec._init,spec._extended) ;
2254 spec._init = kTRUE ;
2255 return ret ;
2256}
2257
2258
2259
2260
2261
2262////////////////////////////////////////////////////////////////////////////////
2263/// Generate a new dataset containing the specified variables with
2264/// events sampled from our distribution.
2265///
2266/// \param[in] whatVars Generate a dataset with the variables (and categories) in this set.
2267/// Any variables of this PDF that are not in `whatVars` will use their
2268/// current values and be treated as fixed parameters.
2269/// \param[in] nEvents Generate the specified number of events or else try to use
2270/// expectedEvents() if nEvents <= 0 (default).
2271/// \param[in] verbose Show which generator strategies are being used.
2272/// \param[in] autoBinned If original distribution is binned, return bin centers and randomise weights
2273/// instead of generating single events.
2274/// \param[in] binnedTag
2275/// \param[in] expectedData Call setExpectedData on the genContext.
2276/// \param[in] extended Randomise number of events generated according to Poisson(nEvents). Only useful
2277/// if PDF is extended.
2278/// \return New dataset. Returns zero in case of an error. The caller takes ownership of the returned
2279/// dataset.
2280
2281RooDataSet *RooAbsPdf::generate(const RooArgSet &whatVars, Double_t nEvents, Bool_t verbose, Bool_t autoBinned, const char* binnedTag, Bool_t expectedData, Bool_t extended) const
2282{
2283 if (nEvents==0 && extendMode()==CanNotBeExtended) {
2284 return new RooDataSet("emptyData","emptyData",whatVars) ;
2285 }
2286
2287 // Request for binned generation
2288 RooAbsGenContext *context = autoGenContext(whatVars,0,0,verbose,autoBinned,binnedTag) ;
2289 if (expectedData) {
2290 context->setExpectedData(kTRUE) ;
2291 }
2292
2293 RooDataSet *generated = 0;
2294 if(0 != context && context->isValid()) {
2295 generated= context->generate(nEvents, kFALSE, extended);
2296 }
2297 else {
2298 coutE(Generation) << "RooAbsPdf::generate(" << GetName() << ") cannot create a valid context" << endl;
2299 }
2300 if(0 != context) delete context;
2301 return generated;
2302}
2303
2304
2305
2306
2307////////////////////////////////////////////////////////////////////////////////
2308/// Internal method
2309
2310RooDataSet *RooAbsPdf::generate(RooAbsGenContext& context, const RooArgSet &whatVars, const RooDataSet *prototype,
2311 Double_t nEvents, Bool_t /*verbose*/, Bool_t randProtoOrder, Bool_t resampleProto,
2312 Bool_t skipInit, Bool_t extended) const
2313{
2314 if (nEvents==0 && (prototype==0 || prototype->numEntries()==0)) {
2315 return new RooDataSet("emptyData","emptyData",whatVars) ;
2316 }
2317
2318 RooDataSet *generated = 0;
2319
2320 // Resampling implies reshuffling in the implementation
2321 if (resampleProto) {
2322 randProtoOrder=kTRUE ;
2323 }
2324
2325 if (randProtoOrder && prototype && prototype->numEntries()!=nEvents) {
2326 coutI(Generation) << "RooAbsPdf::generate (Re)randomizing event order in prototype dataset (Nevt=" << nEvents << ")" << endl ;
2327 Int_t* newOrder = randomizeProtoOrder(prototype->numEntries(),Int_t(nEvents),resampleProto) ;
2328 context.setProtoDataOrder(newOrder) ;
2329 delete[] newOrder ;
2330 }
2331
2332 if(context.isValid()) {
2333 generated= context.generate(nEvents,skipInit,extended);
2334 }
2335 else {
2336 coutE(Generation) << "RooAbsPdf::generate(" << GetName() << ") do not have a valid generator context" << endl;
2337 }
2338 return generated;
2339}
2340
2341
2342
2343
2344////////////////////////////////////////////////////////////////////////////////
2345/// Generate a new dataset using a prototype dataset as a model,
2346/// with values of the variables in `whatVars` sampled from our distribution.
2347///
2348/// \param[in] whatVars Generate for these variables.
2349/// \param[in] prototype Use this dataset
2350/// as a prototype: the new dataset will contain the same number of
2351/// events as the prototype (by default), and any prototype variables not in
2352/// whatVars will be copied into the new dataset for each generated
2353/// event and also used to set our PDF parameters. The user can specify a
2354/// number of events to generate that will override the default. The result is a
2355/// copy of the prototype dataset with only variables in whatVars
2356/// randomized. Variables in whatVars that are not in the prototype
2357/// will be added as new columns to the generated dataset.
2358/// \param[in] nEvents Number of events to generate. Defaults to 0, which means number
2359/// of event in prototype dataset.
2360/// \param[in] verbose Show which generator strategies are being used.
2361/// \param[in] randProtoOrder Randomise order of retrieval of events from proto dataset.
2362/// \param[in] resampleProto Resample from the proto dataset.
2363/// \return The new dataset. Returns zero in case of an error. The caller takes ownership of the
2364/// returned dataset.
2365
2366RooDataSet *RooAbsPdf::generate(const RooArgSet &whatVars, const RooDataSet& prototype,
2367 Int_t nEvents, Bool_t verbose, Bool_t randProtoOrder, Bool_t resampleProto) const
2368{
2369 RooAbsGenContext *context= genContext(whatVars,&prototype,0,verbose);
2370 if (context) {
2371 RooDataSet* data = generate(*context,whatVars,&prototype,nEvents,verbose,randProtoOrder,resampleProto) ;
2372 delete context ;
2373 return data ;
2374 } else {
2375 coutE(Generation) << "RooAbsPdf::generate(" << GetName() << ") ERROR creating generator context" << endl ;
2376 return 0 ;
2377 }
2378}
2379
2380
2381
2382////////////////////////////////////////////////////////////////////////////////
2383/// Return lookup table with randomized access order for prototype events,
2384/// given nProto prototype data events and nGen events that will actually
2385/// be accessed
2386
2388{
2389 // Make unsorted linked list of indeces
2391 Int_t i ;
2392 for (i=0 ; i<nProto ; i++) {
2393 l.Add(new RooInt(i)) ;
2394 }
2395
2396 // Make output list
2397 Int_t* lut = new Int_t[nProto] ;
2398
2399 // Randomly samply input list into output list
2400 if (!resampleProto) {
2401 // In this mode, randomization is a strict reshuffle of the order
2402 for (i=0 ; i<nProto ; i++) {
2403 Int_t iran = RooRandom::integer(nProto-i) ;
2404 RooInt* sample = (RooInt*) l.At(iran) ;
2405 lut[i] = *sample ;
2406 l.Remove(sample) ;
2407 delete sample ;
2408 }
2409 } else {
2410 // In this mode, we resample, i.e. events can be used more than once
2411 for (i=0 ; i<nProto ; i++) {
2412 lut[i] = RooRandom::integer(nProto);
2413 }
2414 }
2415
2416
2417 return lut ;
2418}
2419
2420
2421
2422////////////////////////////////////////////////////////////////////////////////
2423/// Load generatedVars with the subset of directVars that we can generate events for,
2424/// and return a code that specifies the generator algorithm we will use. A code of
2425/// zero indicates that we cannot generate any of the directVars (in this case, nothing
2426/// should be added to generatedVars). Any non-zero codes will be passed to our generateEvent()
2427/// implementation, but otherwise its value is arbitrary. The default implemetation of
2428/// this method returns zero. Subclasses will usually implement this method using the
2429/// matchArgs() methods to advertise the algorithms they provide.
2430
2431Int_t RooAbsPdf::getGenerator(const RooArgSet &/*directVars*/, RooArgSet &/*generatedVars*/, Bool_t /*staticInitOK*/) const
2432{
2433 return 0 ;
2434}
2435
2436
2437
2438////////////////////////////////////////////////////////////////////////////////
2439/// Interface for one-time initialization to setup the generator for the specified code.
2440
2442{
2443}
2444
2445
2446
2447////////////////////////////////////////////////////////////////////////////////
2448/// Interface for generation of an event using the algorithm
2449/// corresponding to the specified code. The meaning of each code is
2450/// defined by the getGenerator() implementation. The default
2451/// implementation does nothing.
2452
2454{
2455}
2456
2457
2458
2459////////////////////////////////////////////////////////////////////////////////
2460/// Check if given observable can be safely generated using the
2461/// pdfs internal generator mechanism (if that existsP). Observables
2462/// on which a PDF depends via more than route are not safe
2463/// for use with internal generators because they introduce
2464/// correlations not known to the internal generator
2465
2467{
2468 // Arg must be direct server of self
2469 if (!findServer(arg.GetName())) return kFALSE ;
2470
2471 // There must be no other dependency routes
2472 for (const auto server : _serverList) {
2473 if(server == &arg) continue;
2474 if(server->dependsOn(arg)) {
2475 return kFALSE ;
2476 }
2477 }
2478
2479 return kTRUE ;
2480}
2481
2482
2483////////////////////////////////////////////////////////////////////////////////
2484/// Generate a new dataset containing the specified variables with events sampled from our distribution.
2485/// \param[in] whatVars Choose variables in which to generate events. Variables not listed here will remain
2486/// constant and not be used for event generation
2487/// \param[in] arg1 Optional RooCmdArg to change behaviour of generateBinned()
2488/// \return RooDataHist *, to be managed by caller.
2489///
2490/// Generate the specified number of events or expectedEvents() if not specified.
2491///
2492/// Any variables of this PDF that are not in whatVars will use their
2493/// current values and be treated as fixed parameters. Returns zero
2494/// in case of an error. The caller takes ownership of the returned
2495/// dataset.
2496///
2497/// The following named arguments are supported
2498/// | Type of CmdArg | Effect on generation
2499/// |-------------------------|-----------------------
2500/// | `Name(const char* name)` | Name of the output dataset
2501/// | `Verbose(Bool_t flag)` | Print informational messages during event generation
2502/// | `NumEvent(int nevt)` | Generate specified number of events
2503/// | `Extended()` | The actual number of events generated will be sampled from a Poisson distribution with mu=nevt. For use with extended maximum likelihood fits
2504/// | `ExpectedData()` | Return a binned dataset _without_ statistical fluctuations (also aliased as Asimov())
2505///
2506
2507RooDataHist *RooAbsPdf::generateBinned(const RooArgSet& whatVars, const RooCmdArg& arg1,const RooCmdArg& arg2,
2508 const RooCmdArg& arg3,const RooCmdArg& arg4, const RooCmdArg& arg5,const RooCmdArg& arg6) const
2509{
2510
2511 // Select the pdf-specific commands
2512 RooCmdConfig pc(Form("RooAbsPdf::generate(%s)",GetName())) ;
2513 pc.defineString("dsetName","Name",0,"") ;
2514 pc.defineInt("verbose","Verbose",0,0) ;
2515 pc.defineInt("extended","Extended",0,0) ;
2516 pc.defineInt("nEvents","NumEvents",0,0) ;
2517 pc.defineDouble("nEventsD","NumEventsD",0,-1.) ;
2518 pc.defineInt("expectedData","ExpectedData",0,0) ;
2519
2520 // Process and check varargs
2521 pc.process(arg1,arg2,arg3,arg4,arg5,arg6) ;
2522 if (!pc.ok(kTRUE)) {
2523 return 0 ;
2524 }
2525
2526 // Decode command line arguments
2527 Double_t nEvents = pc.getDouble("nEventsD") ;
2528 if (nEvents<0) {
2529 nEvents = pc.getInt("nEvents") ;
2530 }
2531 //Bool_t verbose = pc.getInt("verbose") ;
2532 Bool_t extended = pc.getInt("extended") ;
2533 Bool_t expectedData = pc.getInt("expectedData") ;
2534 const char* dsetName = pc.getString("dsetName") ;
2535
2536 if (extended) {
2537 //nEvents = (nEvents==0?Int_t(expectedEvents(&whatVars)+0.5):nEvents) ;
2538 nEvents = (nEvents==0 ? expectedEvents(&whatVars) :nEvents) ;
2539 cxcoutI(Generation) << " Extended mode active, number of events generated (" << nEvents << ") is Poisson fluctuation on "
2540 << GetName() << "::expectedEvents() = " << nEvents << endl ;
2541 // If Poisson fluctuation results in zero events, stop here
2542 if (nEvents==0) {
2543 return 0 ;
2544 }
2545 } else if (nEvents==0) {
2546 cxcoutI(Generation) << "No number of events specified , number of events generated is "
2547 << GetName() << "::expectedEvents() = " << expectedEvents(&whatVars)<< endl ;
2548 }
2549
2550 // Forward to appropriate implementation
2551 RooDataHist* data = generateBinned(whatVars,nEvents,expectedData,extended) ;
2552
2553 // Rename dataset to given name if supplied
2554 if (dsetName && strlen(dsetName)>0) {
2555 data->SetName(dsetName) ;
2556 }
2557
2558 return data ;
2559}
2560
2561
2562
2563
2564////////////////////////////////////////////////////////////////////////////////
2565/// Generate a new dataset containing the specified variables with
2566/// events sampled from our distribution. Generate the specified
2567/// number of events or else try to use expectedEvents() if nEvents <= 0.
2568///
2569/// If expectedData is kTRUE (it is kFALSE by default), the returned histogram returns the 'expected'
2570/// data sample, i.e. no statistical fluctuations are present.
2571///
2572/// Any variables of this PDF that are not in whatVars will use their
2573/// current values and be treated as fixed parameters. Returns zero
2574/// in case of an error. The caller takes ownership of the returned
2575/// dataset.
2576
2577RooDataHist *RooAbsPdf::generateBinned(const RooArgSet &whatVars, Double_t nEvents, Bool_t expectedData, Bool_t extended) const
2578{
2579 // Create empty RooDataHist
2580 RooDataHist* hist = new RooDataHist("genData","genData",whatVars) ;
2581
2582 // Scale to number of events and introduce Poisson fluctuations
2583 if (nEvents<=0) {
2584 if (!canBeExtended()) {
2585 coutE(InputArguments) << "RooAbsPdf::generateBinned(" << GetName() << ") ERROR: No event count provided and p.d.f does not provide expected number of events" << endl ;
2586 delete hist ;
2587 return 0 ;
2588 } else {
2589
2590 // Don't round in expectedData or extended mode
2591 if (expectedData || extended) {
2592 nEvents = expectedEvents(&whatVars) ;
2593 } else {
2594 nEvents = Int_t(expectedEvents(&whatVars)+0.5) ;
2595 }
2596 }
2597 }
2598
2599 // Sample p.d.f. distribution
2600 fillDataHist(hist,&whatVars,1,kTRUE) ;
2601
2602 vector<int> histOut(hist->numEntries()) ;
2603 Double_t histMax(-1) ;
2604 Int_t histOutSum(0) ;
2605 for (int i=0 ; i<hist->numEntries() ; i++) {
2606 hist->get(i) ;
2607 if (expectedData) {
2608
2609 // Expected data, multiply p.d.f by nEvents
2610 Double_t w=hist->weight()*nEvents ;
2611 hist->set(w,sqrt(w)) ;
2612
2613 } else if (extended) {
2614
2615 // Extended mode, set contents to Poisson(pdf*nEvents)
2616 Double_t w = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ;
2617 hist->set(w,sqrt(w)) ;
2618
2619 } else {
2620
2621 // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill
2622 // histogram yet.
2623 if (hist->weight()>histMax) {
2624 histMax = hist->weight() ;
2625 }
2626 histOut[i] = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ;
2627 histOutSum += histOut[i] ;
2628 }
2629 }
2630
2631
2632 if (!expectedData && !extended) {
2633
2634 // Second pass for regular mode - Trim/Extend dataset to exact number of entries
2635
2636 // Calculate difference between what is generated so far and what is requested
2637 Int_t nEvtExtra = abs(Int_t(nEvents)-histOutSum) ;
2638 Int_t wgt = (histOutSum>nEvents) ? -1 : 1 ;
2639
2640 // Perform simple binned accept/reject procedure to get to exact event count
2641 while(nEvtExtra>0) {
2642
2643 Int_t ibinRand = RooRandom::randomGenerator()->Integer(hist->numEntries()) ;
2644 hist->get(ibinRand) ;
2645 Double_t ranY = RooRandom::randomGenerator()->Uniform(histMax) ;
2646
2647 if (ranY<hist->weight()) {
2648 if (wgt==1) {
2649 histOut[ibinRand]++ ;
2650 } else {
2651 // If weight is negative, prior bin content must be at least 1
2652 if (histOut[ibinRand]>0) {
2653 histOut[ibinRand]-- ;
2654 } else {
2655 continue ;
2656 }
2657 }
2658 nEvtExtra-- ;
2659 }
2660 }
2661
2662 // Transfer working array to histogram
2663 for (int i=0 ; i<hist->numEntries() ; i++) {
2664 hist->get(i) ;
2665 hist->set(histOut[i],sqrt(1.0*histOut[i])) ;
2666 }
2667
2668 } else if (expectedData) {
2669
2670 // Second pass for expectedData mode -- Normalize to exact number of requested events
2671 // Minor difference may be present in first round due to difference between
2672 // bin average and bin integral in sampling bins
2673 Double_t corr = nEvents/hist->sumEntries() ;
2674 for (int i=0 ; i<hist->numEntries() ; i++) {
2675 hist->get(i) ;
2676 hist->set(hist->weight()*corr,sqrt(hist->weight()*corr)) ;
2677 }
2678
2679 }
2680
2681 return hist;
2682}
2683
2684
2685
2686////////////////////////////////////////////////////////////////////////////////
2687/// Special generator interface for generation of 'global observables' -- for RooStats tools
2688
2690{
2691 return generate(whatVars,nEvents) ;
2692}
2693
2694namespace {
2695void removeRangeOverlap(std::vector<std::pair<double, double>>& ranges) {
2696 //Sort from left to right
2697 std::sort(ranges.begin(), ranges.end());
2698
2699 for (auto it = ranges.begin(); it != ranges.end(); ++it) {
2700 double& startL = it->first;
2701 double& endL = it->second;
2702
2703 for (auto innerIt = it+1; innerIt != ranges.end(); ++innerIt) {
2704 const double startR = innerIt->first;
2705 const double endR = innerIt->second;
2706
2707 if (startL <= startR && startR <= endL) {
2708 //Overlapping ranges, extend left one
2709 endL = std::max(endL, endR);
2710 *innerIt = make_pair(0., 0.);
2711 }
2712 }
2713 }
2714
2715 auto newEnd = std::remove_if(ranges.begin(), ranges.end(),
2716 [](const std::pair<double,double>& input){
2717 return input.first == input.second;
2718 });
2719 ranges.erase(newEnd, ranges.end());
2720}
2721}
2722
2723
2724////////////////////////////////////////////////////////////////////////////////
2725/// Plot (project) PDF on specified frame. If a PDF is plotted in an empty frame, it
2726/// will show a unit normalized curve in the frame variable, taken at the present value
2727/// of other observables defined for this PDF.
2728///
2729/// If a PDF is plotted in a frame in which a dataset has already been plotted, it will
2730/// show a projected curve integrated over all variables that were present in the shown
2731/// dataset except for the one on the x-axis. The normalization of the curve will also
2732/// be adjusted to the event count of the plotted dataset. An informational message
2733/// will be printed for each projection step that is performed.
2734///
2735/// This function takes the following named arguments (for more arguments, see also
2736/// RooAbsReal::plotOn(RooPlot*,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,
2737/// const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,
2738/// const RooCmdArg&) const )
2739///
2740///
2741/// <table>
2742/// <tr><th> Type of CmdArg <th> Projection control
2743/// <tr><td> `Slice(const RooArgSet& set)` <td> Override default projection behaviour by omitting
2744/// observables listed in set from the projection, resulting a 'slice' plot. Slicing is usually
2745/// only sensible in discrete observables
2746/// <tr><td> `Project(const RooArgSet& set)` <td> Override default projection behaviour by projecting
2747/// over observables given in set and complete ignoring the default projection behavior. Advanced use only.
2748/// <tr><td> `ProjWData(const RooAbsData& d)` <td> Override default projection _technique_ (integration). For observables
2749/// present in given dataset projection of PDF is achieved by constructing an average over all observable
2750/// values in given set. Consult RooFit plotting tutorial for further explanation of meaning & use of this technique
2751/// <tr><td> `ProjWData(const RooArgSet& s, const RooAbsData& d)` <td> As above but only consider subset 's' of
2752/// observables in dataset 'd' for projection through data averaging
2753/// <tr><td> `ProjectionRange(const char* rn)` <td> Override default range of projection integrals to a different
2754/// range specified by given range name. This technique allows you to project a finite width slice in a real-valued observable
2755/// <tr><td> `NormRange(const char* name)` <td> Calculate curve normalization w.r.t. specified range[s].
2756/// \note A Range() by default sets a NormRange() on the same range, but this option allows
2757/// to override the default, or specify normalization ranges when the full curve is to be drawn.
2758///
2759///
2760/// <tr><th><th> Misc content control
2761/// <tr><td> `Normalization(Double_t scale, ScaleType code)` <td> Adjust normalization by given scale factor.
2762/// Interpretation of number depends on code:
2763/// `RooAbsReal::Relative`: relative adjustment factor
2764/// `RooAbsReal::NumEvent`: scale to match given number of events.
2765/// <tr><td> `Name(const chat* name)` <td> Give curve specified name in frame. Useful if curve is to be referenced later
2766/// <tr><td> `Asymmetry(const RooCategory& c)` <td> Show the asymmetry of the PDF in given two-state category
2767/// \f$ \frac{F(+)-F(-)}{F(+)+F(-)} \f$ rather than the PDF projection. Category must have two
2768/// states with indices -1 and +1 or three states with indeces -1,0 and +1.
2769/// <tr><td> `ShiftToZero(Bool_t flag)` <td> Shift entire curve such that lowest visible point is at exactly zero.
2770/// Mostly useful when plotting -log(L) or \f$ \chi^2 \f$ distributions
2771/// <tr><td> `AddTo(const char* name, double_t wgtSelf, double_t wgtOther)` <td> Add constructed projection to
2772/// already existing curve with given name and relative weight factors
2773/// <tr><td> `Components(const char* names)` <td> When plotting sums of PDFs, plot only the named components (*e.g.* only
2774/// the signal of a signal+background model).
2775/// <tr><td> `Components(const RooArgSet& compSet)` <td> As above, but pass a RooArgSet of the components themselves.
2776///
2777/// <tr><th><th> Plotting control
2778/// <tr><td> `LineStyle(Int_t style)` <td> Select line style by ROOT line style code, default is solid
2779/// <tr><td> `LineColor(Int_t color)` <td> Select line color by ROOT color code, default is blue
2780/// <tr><td> `LineWidth(Int_t width)` <td> Select line with in pixels, default is 3
2781/// <tr><td> `FillStyle(Int_t style)` <td> Select fill style, default is not filled. If a filled style is selected,
2782/// also use VLines() to add vertical downward lines at end of curve to ensure proper closure
2783/// <tr><td> `FillColor(Int_t color)` <td> Select fill color by ROOT color code
2784/// <tr><td> `Range(const char* name)` <td> Only draw curve in range defined by given name. Multiple comma-separated ranges can be given.
2785/// <tr><td> `Range(double lo, double hi)` <td> Only draw curve in specified range
2786/// <tr><td> `VLines()` <td> Add vertical lines to y=0 at end points of curve
2787/// <tr><td> `Precision(Double_t eps)` <td> Control precision of drawn curve w.r.t to scale of plot, default is 1e-3. A higher precision will
2788/// result in more and more densely spaced curve points. A negative precision value will disable
2789/// adaptive point spacing and restrict sampling to the grid point of points defined by the binning
2790/// of the plotted observable (recommended for expensive functions such as profile likelihoods)
2791/// <tr><td> `Invisible(Bool_t flag)` <td> Add curve to frame, but do not display. Useful in combination AddTo()
2792/// <tr><td> `VisualizeError(const RooFitResult& fitres, Double_t Z=1, Bool_t linearMethod=kTRUE)`
2793/// <td> Visualize the uncertainty on the parameters, as given in fitres, at 'Z' sigma.
2794///
2795/// <tr><td> `VisualizeError(const RooFitResult& fitres, const RooArgSet& param, Double_t Z=1, Bool_t linearMethod=kTRUE)`
2796/// <td> Visualize the uncertainty on the subset of parameters 'param', as given in fitres, at 'Z' sigma.
2797/// </table>
2798
2800{
2801
2802 // Pre-processing if p.d.f. contains a fit range and there is no command specifying one,
2803 // add a fit range as default range
2804 RooCmdArg* plotRange(0) ;
2805 RooCmdArg* normRange2(0) ;
2806 if (getStringAttribute("fitrange") && !cmdList.FindObject("Range") &&
2807 !cmdList.FindObject("RangeWithName")) {
2808 plotRange = (RooCmdArg*) RooFit::Range(getStringAttribute("fitrange")).Clone() ;
2809 cmdList.Add(plotRange) ;
2810 }
2811
2812 if (getStringAttribute("fitrange") && !cmdList.FindObject("NormRange")) {
2813 normRange2 = (RooCmdArg*) RooFit::NormRange(getStringAttribute("fitrange")).Clone() ;
2814 cmdList.Add(normRange2) ;
2815 }
2816
2817 if (plotRange || normRange2) {
2818 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") p.d.f was fitted in range and no explicit "
2819 << (plotRange?"plot":"") << ((plotRange&&normRange2)?",":"")
2820 << (normRange2?"norm":"") << " range was specified, using fit range as default" << endl ;
2821 }
2822
2823 // Sanity checks
2824 if (plotSanityChecks(frame)) return frame ;
2825
2826 // Select the pdf-specific commands
2827 RooCmdConfig pc(Form("RooAbsPdf::plotOn(%s)",GetName())) ;
2828 pc.defineDouble("scaleFactor","Normalization",0,1.0) ;
2829 pc.defineInt("scaleType","Normalization",0,Relative) ;
2830 pc.defineObject("compSet","SelectCompSet",0) ;
2831 pc.defineString("compSpec","SelectCompSpec",0) ;
2832 pc.defineObject("asymCat","Asymmetry",0) ;
2833 pc.defineDouble("rangeLo","Range",0,-999.) ;
2834 pc.defineDouble("rangeHi","Range",1,-999.) ;
2835 pc.defineString("rangeName","RangeWithName",0,"") ;
2836 pc.defineString("normRangeName","NormRange",0,"") ;
2837 pc.defineInt("rangeAdjustNorm","Range",0,0) ;
2838 pc.defineInt("rangeWNAdjustNorm","RangeWithName",0,0) ;
2839 pc.defineMutex("SelectCompSet","SelectCompSpec") ;
2840 pc.defineMutex("Range","RangeWithName") ;
2841 pc.allowUndefined() ; // unknowns may be handled by RooAbsReal
2842
2843 // Process and check varargs
2844 pc.process(cmdList) ;
2845 if (!pc.ok(kTRUE)) {
2846 return frame ;
2847 }
2848
2849 // Decode command line arguments
2850 ScaleType stype = (ScaleType) pc.getInt("scaleType") ;
2851 Double_t scaleFactor = pc.getDouble("scaleFactor") ;
2852 const RooAbsCategoryLValue* asymCat = (const RooAbsCategoryLValue*) pc.getObject("asymCat") ;
2853 const char* compSpec = pc.getString("compSpec") ;
2854 const RooArgSet* compSet = (const RooArgSet*) pc.getObject("compSet") ;
2855 Bool_t haveCompSel = ( (compSpec && strlen(compSpec)>0) || compSet) ;
2856
2857 // Suffix for curve name
2858 TString nameSuffix ;
2859 if (compSpec && strlen(compSpec)>0) {
2860 nameSuffix.Append("_Comp[") ;
2861 nameSuffix.Append(compSpec) ;
2862 nameSuffix.Append("]") ;
2863 } else if (compSet) {
2864 nameSuffix.Append("_Comp[") ;
2865 nameSuffix.Append(compSet->contentsString().c_str()) ;
2866 nameSuffix.Append("]") ;
2867 }
2868
2869 // Remove PDF-only commands from command list
2870 pc.stripCmdList(cmdList,"SelectCompSet,SelectCompSpec") ;
2871
2872 // Adjust normalization, if so requested
2873 if (asymCat) {
2874 RooCmdArg cnsuffix("CurveNameSuffix",0,0,0,0,nameSuffix.Data(),0,0,0) ;
2875 cmdList.Add(&cnsuffix);
2876 return RooAbsReal::plotOn(frame,cmdList) ;
2877 }
2878
2879 // More sanity checks
2880 Double_t nExpected(1) ;
2881 if (stype==RelativeExpected) {
2882 if (!canBeExtended()) {
2883 coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName()
2884 << "): ERROR the 'Expected' scale option can only be used on extendable PDFs" << endl ;
2885 return frame ;
2886 }
2887 nExpected = expectedEvents(frame->getNormVars()) ;
2888 }
2889
2890 if (stype != Raw) {
2891
2892 if (frame->getFitRangeNEvt() && stype==Relative) {
2893
2894 Bool_t hasCustomRange(kFALSE), adjustNorm(kFALSE) ;
2895
2896 std::vector<pair<Double_t,Double_t> > rangeLim;
2897
2898 // Retrieve plot range to be able to adjust normalization to data
2899 if (pc.hasProcessed("Range")) {
2900
2901 Double_t rangeLo = pc.getDouble("rangeLo") ;
2902 Double_t rangeHi = pc.getDouble("rangeHi") ;
2903 rangeLim.push_back(make_pair(rangeLo,rangeHi)) ;
2904 adjustNorm = pc.getInt("rangeAdjustNorm") ;
2905 hasCustomRange = kTRUE ;
2906
2907 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") only plotting range ["
2908 << rangeLo << "," << rangeHi << "]" ;
2909 if (!pc.hasProcessed("NormRange")) {
2910 ccoutI(Plotting) << ", curve is normalized to data in " << (adjustNorm?"given":"full") << " range" << endl ;
2911 } else {
2912 ccoutI(Plotting) << endl ;
2913 }
2914
2915 nameSuffix.Append(Form("_Range[%f_%f]",rangeLo,rangeHi)) ;
2916
2917 } else if (pc.hasProcessed("RangeWithName")) {
2918 for (const std::string& rangeNameToken : RooHelpers::tokenise(pc.getString("rangeName",0,true), ",")) {
2919 if (!frame->getPlotVar()->hasRange(rangeNameToken.c_str())) {
2920 coutE(Plotting) << "Range '" << rangeNameToken << "' not defined for variable '"
2921 << frame->getPlotVar()->GetName() << "'. Ignoring ..." << std::endl;
2922 continue;
2923 }
2924 const double rangeLo = frame->getPlotVar()->getMin(rangeNameToken.c_str());
2925 const double rangeHi = frame->getPlotVar()->getMax(rangeNameToken.c_str());
2926 rangeLim.push_back(make_pair(rangeLo,rangeHi));
2927 }
2928 adjustNorm = pc.getInt("rangeWNAdjustNorm") ;
2929 hasCustomRange = kTRUE ;
2930
2931 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") only plotting range '" << pc.getString("rangeName",0,kTRUE) << "'" ;
2932 if (!pc.hasProcessed("NormRange")) {
2933 ccoutI(Plotting) << ", curve is normalized to data in " << (adjustNorm?"given":"full") << " range" << endl ;
2934 } else {
2935 ccoutI(Plotting) << endl ;
2936 }
2937
2938 nameSuffix.Append(Form("_Range[%s]",pc.getString("rangeName"))) ;
2939 }
2940 // Specification of a normalization range override those in a regular range
2941 if (pc.hasProcessed("NormRange")) {
2942 rangeLim.clear();
2943 for (const auto& rangeNameToken : RooHelpers::tokenise(pc.getString("normRangeName",0,true), ",")) {
2944 if (!frame->getPlotVar()->hasRange(rangeNameToken.c_str())) {
2945 coutE(Plotting) << "Range '" << rangeNameToken << "' not defined for variable '"
2946 << frame->getPlotVar()->GetName() << "'. Ignoring ..." << std::endl;
2947 continue;
2948 }
2949 const double rangeLo = frame->getPlotVar()->getMin(rangeNameToken.c_str());
2950 const double rangeHi = frame->getPlotVar()->getMax(rangeNameToken.c_str());
2951 rangeLim.push_back(make_pair(rangeLo,rangeHi));
2952 }
2953 adjustNorm = kTRUE ;
2954 hasCustomRange = kTRUE ;
2955 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") p.d.f. curve is normalized using explicit choice of ranges '" << pc.getString("normRangeName",0,kTRUE) << "'" << endl ;
2956
2957 nameSuffix.Append(Form("_NormRange[%s]",pc.getString("rangeName"))) ;
2958
2959 }
2960
2961 if (hasCustomRange && adjustNorm) {
2962 const std::size_t oldSize = rangeLim.size();
2963 removeRangeOverlap(rangeLim);
2964
2965 if (oldSize != rangeLim.size()) {
2966 // User gave overlapping ranges. This leads to double-counting events and integrals, and must
2967 // therefore be avoided.
2968 coutE(Plotting) << "Requested ranges overlap. For correct plotting, new ranges "
2969 "will be defined." << std::endl;
2970 auto plotVar = dynamic_cast<RooRealVar*>(frame->getPlotVar());
2971 assert(plotVar);
2972 std::string rangesNoOverlap;
2973 for (auto it = rangeLim.begin(); it != rangeLim.end(); ++it) {
2974 std::stringstream rangeName;
2975 rangeName << "Remove_overlap_range_" << it - rangeLim.begin();
2976 plotVar->setRange(rangeName.str().c_str(), it->first, it->second);
2977 if (!rangesNoOverlap.empty())
2978 rangesNoOverlap += ",";
2979 rangesNoOverlap += rangeName.str();
2980 }
2981
2982 auto rangeArg = static_cast<RooCmdArg*>(cmdList.FindObject("RangeWithName"));
2983 if (rangeArg)
2984 rangeArg->setString(0, rangesNoOverlap.c_str());
2985 else {
2986 plotRange = new RooCmdArg(RooFit::Range(rangesNoOverlap.c_str()));
2987 cmdList.Add(plotRange);
2988 }
2989 }
2990
2991 Double_t rangeNevt(0) ;
2992 for (const auto& riter : rangeLim) {
2993 Double_t nevt= frame->getFitRangeNEvt(riter.first, riter.second);
2994 rangeNevt += nevt ;
2995 }
2996
2997 scaleFactor *= rangeNevt/nExpected ;
2998
2999 } else {
3000 scaleFactor *= frame->getFitRangeNEvt()/nExpected ;
3001 }
3002 } else if (stype==RelativeExpected) {
3003 scaleFactor *= nExpected ;
3004 } else if (stype==NumEvent) {
3005 scaleFactor /= nExpected ;
3006 }
3007 scaleFactor *= frame->getFitRangeBinW() ;
3008 }
3009 frame->updateNormVars(*frame->getPlotVar()) ;
3010
3011 // Append overriding scale factor command at end of original command list
3012 RooCmdArg tmp = RooFit::Normalization(scaleFactor,Raw) ;
3013 tmp.setInt(1,1) ; // Flag this normalization command as created for internal use (so that VisualizeError can strip it)
3014 cmdList.Add(&tmp) ;
3015
3016 // Was a component selected requested
3017 if (haveCompSel) {
3018
3019 // Get complete set of tree branch nodes
3020 RooArgSet branchNodeSet ;
3021 branchNodeServerList(&branchNodeSet) ;
3022
3023 // Discard any non-RooAbsReal nodes
3024 for (const auto arg : branchNodeSet) {
3025 if (!dynamic_cast<RooAbsReal*>(arg)) {
3026 branchNodeSet.remove(*arg) ;
3027 }
3028 }
3029
3030 // Obtain direct selection
3031 RooArgSet* dirSelNodes ;
3032 if (compSet) {
3033 dirSelNodes = (RooArgSet*) branchNodeSet.selectCommon(*compSet) ;
3034 } else {
3035 dirSelNodes = (RooArgSet*) branchNodeSet.selectByName(compSpec) ;
3036 }
3037 if (dirSelNodes->getSize()>0) {
3038 coutI(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") directly selected PDF components: " << *dirSelNodes << endl ;
3039
3040 // Do indirect selection and activate both
3041 plotOnCompSelect(dirSelNodes) ;
3042 } else {
3043 if (compSet) {
3044 coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") ERROR: component selection set " << *compSet << " does not match any components of p.d.f." << endl ;
3045 } else {
3046 coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName() << ") ERROR: component selection expression '" << compSpec << "' does not select any components of p.d.f." << endl ;
3047 }
3048 return 0 ;
3049 }
3050
3051 delete dirSelNodes ;
3052 }
3053
3054
3055 RooCmdArg cnsuffix("CurveNameSuffix",0,0,0,0,nameSuffix.Data(),0,0,0) ;
3056 cmdList.Add(&cnsuffix);
3057
3058 RooPlot* ret = RooAbsReal::plotOn(frame,cmdList) ;
3059
3060 // Restore selection status ;
3061 if (haveCompSel) plotOnCompSelect(0) ;
3062
3063 if (plotRange) {
3064 delete plotRange ;
3065 }
3066 if (normRange2) {
3067 delete normRange2 ;
3068 }
3069
3070 return ret ;
3071}
3072
3073
3074//_____________________________________________________________________________
3075/// Plot oneself on 'frame'. In addition to features detailed in RooAbsReal::plotOn(),
3076/// the scale factor for a PDF can be interpreted in three different ways. The interpretation
3077/// is controlled by ScaleType
3078/// ```
3079/// Relative - Scale factor is applied on top of PDF normalization scale factor
3080/// NumEvent - Scale factor is interpreted as a number of events. The surface area
3081/// under the PDF curve will match that of a histogram containing the specified
3082/// number of event
3083/// Raw - Scale factor is applied to the raw (projected) probability density.
3084/// Not too useful, option provided for completeness.
3085/// ```
3086// coverity[PASS_BY_VALUE]
3088{
3089
3090 // Sanity checks
3091 if (plotSanityChecks(frame)) return frame ;
3092
3093 // More sanity checks
3094 Double_t nExpected(1) ;
3095 if (o.stype==RelativeExpected) {
3096 if (!canBeExtended()) {
3097 coutE(Plotting) << "RooAbsPdf::plotOn(" << GetName()
3098 << "): ERROR the 'Expected' scale option can only be used on extendable PDFs" << endl ;
3099 return frame ;
3100 }
3101 nExpected = expectedEvents(frame->getNormVars()) ;
3102 }
3103
3104 // Adjust normalization, if so requested
3105 if (o.stype != Raw) {
3106
3107 if (frame->getFitRangeNEvt() && o.stype==Relative) {
3108 // If non-default plotting range is specified, adjust number of events in fit range
3109 o.scaleFactor *= frame->getFitRangeNEvt()/nExpected ;
3110 } else if (o.stype==RelativeExpected) {
3111 o.scaleFactor *= nExpected ;
3112 } else if (o.stype==NumEvent) {
3113 o.scaleFactor /= nExpected ;
3114 }
3115 o.scaleFactor *= frame->getFitRangeBinW() ;
3116 }
3117 frame->updateNormVars(*frame->getPlotVar()) ;
3118
3119 return RooAbsReal::plotOn(frame,o) ;
3120}
3121
3122
3123
3124
3125////////////////////////////////////////////////////////////////////////////////
3126/// The following named arguments are supported
3127/// <table>
3128/// <tr><th> Type of CmdArg <th> Effect on parameter box
3129/// <tr><td> `Parameters(const RooArgSet& param)` <td> Only the specified subset of parameters will be shown. By default all non-constant parameters are shown.
3130/// <tr><td> `ShowConstants(Bool_t flag)` <td> Also display constant parameters
3131/// <tr><td> `Format(const char* optStr)` <td> \deprecated Classing parameter formatting options, provided for backward compatibility
3132///
3133/// <tr><td> `Format(const char* what,...)` <td> Parameter formatting options.
3134/// | Parameter | Format
3135/// | ---------------------- | --------------------------
3136/// | `const char* what` | Controls what is shown. "N" adds name, "E" adds error, "A" shows asymmetric error, "U" shows unit, "H" hides the value
3137/// | `FixedPrecision(int n)`| Controls precision, set fixed number of digits
3138/// | `AutoPrecision(int n)` | Controls precision. Number of shown digits is calculated from error + n specified additional digits (1 is sensible default)
3139/// <tr><td> `Label(const chat* label)` <td> Add header label to parameter box
3140/// <tr><td> `Layout(Double_t xmin, Double_t xmax, Double_t ymax)` <td> Specify relative position of left/right side of box and top of box.
3141/// Coordinates are given as position on the pad between 0 and 1.
3142/// The lower end of the box is calculated automatically from the number of lines in the box.
3143/// </table>
3144///
3145///
3146/// Example use:
3147/// ```
3148/// pdf.paramOn(frame, Label("fit result"), Format("NEU",AutoPrecision(1)) ) ;
3149/// ```
3150///
3151
3152RooPlot* RooAbsPdf::paramOn(RooPlot* frame, const RooCmdArg& arg1, const RooCmdArg& arg2,
3153 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
3154 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
3155{
3156 // Stuff all arguments in a list
3157 RooLinkedList cmdList;
3158 cmdList.Add(const_cast<RooCmdArg*>(&arg1)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg2)) ;
3159 cmdList.Add(const_cast<RooCmdArg*>(&arg3)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg4)) ;
3160 cmdList.Add(const_cast<RooCmdArg*>(&arg5)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg6)) ;
3161 cmdList.Add(const_cast<RooCmdArg*>(&arg7)) ; cmdList.Add(const_cast<RooCmdArg*>(&arg8)) ;
3162
3163 // Select the pdf-specific commands
3164 RooCmdConfig pc(Form("RooAbsPdf::paramOn(%s)",GetName())) ;
3165 pc.defineString("label","Label",0,"") ;
3166 pc.defineDouble("xmin","Layout",0,0.50) ;
3167 pc.defineDouble("xmax","Layout",1,0.99) ;
3168 pc.defineInt("ymaxi","Layout",0,Int_t(0.95*10000)) ;
3169 pc.defineInt("showc","ShowConstants",0,0) ;
3170 pc.defineObject("params","Parameters",0,0) ;
3171 pc.defineString("formatStr","Format",0,"NELU") ;
3172 pc.defineInt("sigDigit","Format",0,2) ;
3173 pc.defineInt("dummy","FormatArgs",0,0) ;
3174 pc.defineMutex("Format","FormatArgs") ;
3175
3176 // Process and check varargs
3177 pc.process(cmdList) ;
3178 if (!pc.ok(kTRUE)) {
3179 return frame ;
3180 }
3181
3182 const char* label = pc.getString("label") ;
3183 Double_t xmin = pc.getDouble("xmin") ;
3184 Double_t xmax = pc.getDouble("xmax") ;
3185 Double_t ymax = pc.getInt("ymaxi") / 10000. ;
3186 Int_t showc = pc.getInt("showc") ;
3187
3188
3189 const char* formatStr = pc.getString("formatStr") ;
3190 Int_t sigDigit = pc.getInt("sigDigit") ;
3191
3192 // Decode command line arguments
3193 RooArgSet* params = static_cast<RooArgSet*>(pc.getObject("params")) ;
3194 if (!params) {
3195 params = getParameters(frame->getNormVars()) ;
3196 if (pc.hasProcessed("FormatArgs")) {
3197 const RooCmdArg* formatCmd = static_cast<RooCmdArg*>(cmdList.FindObject("FormatArgs")) ;
3198 paramOn(frame,*params,showc,label,0,0,xmin,xmax,ymax,formatCmd) ;
3199 } else {
3200 paramOn(frame,*params,showc,label,sigDigit,formatStr,xmin,xmax,ymax) ;
3201 }
3202 delete params ;
3203 } else {
3204 RooArgSet* pdfParams = getParameters(frame->getNormVars()) ;
3205 RooArgSet* selParams = static_cast<RooArgSet*>(pdfParams->selectCommon(*params)) ;
3206 if (pc.hasProcessed("FormatArgs")) {
3207 const RooCmdArg* formatCmd = static_cast<RooCmdArg*>(cmdList.FindObject("FormatArgs")) ;
3208 paramOn(frame,*selParams,showc,label,0,0,xmin,xmax,ymax,formatCmd) ;
3209 } else {
3210 paramOn(frame,*selParams,showc,label,sigDigit,formatStr,xmin,xmax,ymax) ;
3211 }
3212 delete selParams ;
3213 delete pdfParams ;
3214 }
3215
3216 return frame ;
3217}
3218
3219
3220
3221
3222////////////////////////////////////////////////////////////////////////////////
3223/// \deprecated Obsolete, provided for backward compatibility. Don't use.
3224
3225RooPlot* RooAbsPdf::paramOn(RooPlot* frame, const RooAbsData* data, const char *label,
3226 Int_t sigDigits, Option_t *options, Double_t xmin,
3228{
3229 RooArgSet* params = getParameters(data) ;
3230 TString opts(options) ;
3231 paramOn(frame,*params,opts.Contains("c"),label,sigDigits,options,xmin,xmax,ymax) ;
3232 delete params ;
3233 return frame ;
3234}
3235
3236
3237
3238////////////////////////////////////////////////////////////////////////////////
3239/// Add a text box with the current parameter values and their errors to the frame.
3240/// Observables of this PDF appearing in the 'data' dataset will be omitted.
3241///
3242/// Optional label will be inserted as first line of the text box. Use 'sigDigits'
3243/// to modify the default number of significant digits printed. The 'xmin,xmax,ymax'
3244/// values specify the initial relative position of the text box in the plot frame
3245
3246RooPlot* RooAbsPdf::paramOn(RooPlot* frame, const RooArgSet& params, Bool_t showConstants, const char *label,
3247 Int_t sigDigits, Option_t *options, Double_t xmin,
3248 Double_t xmax ,Double_t ymax, const RooCmdArg* formatCmd)
3249{
3250
3251 // parse the options
3252 TString opts = options;
3253 opts.ToLower();
3254 Bool_t showLabel= (label != 0 && strlen(label) > 0);
3255
3256 // calculate the box's size, adjusting for constant parameters
3257 TIterator* pIter = params.createIterator() ;
3258
3259 Double_t ymin(ymax), dy(0.06);
3260 RooRealVar *var = 0;
3261 while((var=(RooRealVar*)pIter->Next())) {
3262 if(showConstants || !var->isConstant()) ymin-= dy;
3263 }
3264
3265 if(showLabel) ymin-= dy;
3266
3267 // create the box and set its options
3268 TPaveText *box= new TPaveText(xmin,ymax,xmax,ymin,"BRNDC");
3269 if(!box) return 0;
3270 box->SetName(Form("%s_paramBox",GetName())) ;
3271 box->SetFillColor(0);
3272 box->SetBorderSize(1);
3273 box->SetTextAlign(12);
3274 box->SetTextSize(0.04F);
3275 box->SetFillStyle(1001);
3276 box->SetFillColor(0);
3277 //char buffer[512];
3278 pIter->Reset() ;
3279 while((var=(RooRealVar*)pIter->Next())) {
3280 if(var->isConstant() && !showConstants) continue;
3281
3282 TString *formatted= options ? var->format(sigDigits, options) : var->format(*formatCmd) ;
3283 box->AddText(formatted->Data());
3284 delete formatted;
3285 }
3286 // add the optional label if specified
3287 if(showLabel) box->AddText(label);
3288
3289 // Add box to frame
3290 frame->addObject(box) ;
3291
3292 delete pIter ;
3293 return frame ;
3294}
3295
3296
3297
3298
3299////////////////////////////////////////////////////////////////////////////////
3300/// Return expected number of events from this p.d.f for use in extended
3301/// likelihood calculations. This default implementation returns zero
3302
3304{
3305 return 0 ;
3306}
3307
3308
3309
3310////////////////////////////////////////////////////////////////////////////////
3311/// Change global level of verbosity for p.d.f. evaluations
3312
3314{
3315 _verboseEval = stat ;
3316}
3317
3318
3319
3320////////////////////////////////////////////////////////////////////////////////
3321/// Return global level of verbosity for p.d.f. evaluations
3322
3324{
3325 return _verboseEval ;
3326}
3327
3328
3329
3330////////////////////////////////////////////////////////////////////////////////
3331/// Dummy implementation
3332
3334{
3335}
3336
3337
3338
3339////////////////////////////////////////////////////////////////////////////////
3340/// Destructor of normalization cache element. If this element
3341/// provides the 'current' normalization stored in RooAbsPdf::_norm
3342/// zero _norm pointer here before object pointed to is deleted here
3343
3345{
3346 // Zero _norm pointer in RooAbsPdf if it is points to our cache payload
3347 if (_owner) {
3348 RooAbsPdf* pdfOwner = static_cast<RooAbsPdf*>(_owner) ;
3349 if (pdfOwner->_norm == _norm) {
3350 pdfOwner->_norm = 0 ;
3351 }
3352 }
3353
3354 delete _norm ;
3355}
3356
3357
3358
3359////////////////////////////////////////////////////////////////////////////////
3360/// Return a p.d.f that represent a projection of this p.d.f integrated over given observables
3361
3363{
3364 // Construct name for new object
3365 TString name(GetName()) ;
3366 name.Append("_Proj[") ;
3367 if (iset.getSize()>0) {
3368 TIterator* iter = iset.createIterator() ;
3369 RooAbsArg* arg ;
3370 Bool_t first(kTRUE) ;
3371 while((arg=(RooAbsArg*)iter->Next())) {
3372 if (first) {
3373 first=kFALSE ;
3374 } else {
3375 name.Append(",") ;
3376 }
3377 name.Append(arg->GetName()) ;
3378 }
3379 delete iter ;
3380 }
3381 name.Append("]") ;
3382
3383 // Return projected p.d.f.
3384 return new RooProjectedPdf(name.Data(),name.Data(),*this,iset) ;
3385}
3386
3387
3388
3389////////////////////////////////////////////////////////////////////////////////
3390/// Create a cumulative distribution function of this p.d.f in terms
3391/// of the observables listed in iset. If no nset argument is given
3392/// the c.d.f normalization is constructed over the integrated
3393/// observables, so that its maximum value is precisely 1. It is also
3394/// possible to choose a different normalization for
3395/// multi-dimensional p.d.f.s: eg. for a pdf f(x,y,z) one can
3396/// construct a partial cdf c(x,y) that only when integrated itself
3397/// over z results in a maximum value of 1. To construct such a cdf pass
3398/// z as argument to the optional nset argument
3399
3401{
3402 return createCdf(iset,RooFit::SupNormSet(nset)) ;
3403}
3404
3405
3406
3407////////////////////////////////////////////////////////////////////////////////
3408/// Create an object that represents the integral of the function over one or more observables listed in `iset`.
3409/// The actual integration calculation is only performed when the return object is evaluated. The name
3410/// of the integral object is automatically constructed from the name of the input function, the variables
3411/// it integrates and the range integrates over
3412///
3413/// The following named arguments are accepted
3414/// | Type of CmdArg | Effect on CDF
3415/// | ---------------------|-------------------
3416/// | SupNormSet(const RooArgSet&) | Observables over which should be normalized _in addition_ to the integration observables
3417/// | ScanNumCdf() | Apply scanning technique if cdf integral involves numeric integration [ default ]
3418/// | ScanAllCdf() | Always apply scanning technique
3419/// | ScanNoCdf() | Never apply scanning technique
3420/// | 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
3421
3422RooAbsReal* RooAbsPdf::createCdf(const RooArgSet& iset, const RooCmdArg& arg1, const RooCmdArg& arg2,
3423 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
3424 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8)
3425{
3426 // Define configuration for this method
3427 RooCmdConfig pc(Form("RooAbsReal::createCdf(%s)",GetName())) ;
3428 pc.defineObject("supNormSet","SupNormSet",0,0) ;
3429 pc.defineInt("numScanBins","ScanParameters",0,1000) ;
3430 pc.defineInt("intOrder","ScanParameters",1,2) ;
3431 pc.defineInt("doScanNum","ScanNumCdf",0,1) ;
3432 pc.defineInt("doScanAll","ScanAllCdf",0,0) ;
3433 pc.defineInt("doScanNon","ScanNoCdf",0,0) ;
3434 pc.defineMutex("ScanNumCdf","ScanAllCdf","ScanNoCdf") ;
3435
3436 // Process & check varargs
3437 pc.process(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8) ;
3438 if (!pc.ok(kTRUE)) {
3439 return 0 ;
3440 }
3441
3442 // Extract values from named arguments
3443 const RooArgSet* snset = static_cast<const RooArgSet*>(pc.getObject("supNormSet",0)) ;
3444 RooArgSet nset ;
3445 if (snset) {
3446 nset.add(*snset) ;
3447 }
3448 Int_t numScanBins = pc.getInt("numScanBins") ;
3449 Int_t intOrder = pc.getInt("intOrder") ;
3450 Int_t doScanNum = pc.getInt("doScanNum") ;
3451 Int_t doScanAll = pc.getInt("doScanAll") ;
3452 Int_t doScanNon = pc.getInt("doScanNon") ;
3453
3454 // If scanning technique is not requested make integral-based cdf and return
3455 if (doScanNon) {
3456 return createIntRI(iset,nset) ;
3457 }
3458 if (doScanAll) {
3459 return createScanCdf(iset,nset,numScanBins,intOrder) ;
3460 }
3461 if (doScanNum) {
3463 Int_t isNum= (tmp->numIntRealVars().getSize()>0) ;
3464 delete tmp ;
3465
3466 if (isNum) {
3467 coutI(NumIntegration) << "RooAbsPdf::createCdf(" << GetName() << ") integration over observable(s) " << iset << " involves numeric integration," << endl
3468 << " constructing cdf though numeric integration of sampled pdf in " << numScanBins << " bins and applying order "
3469 << intOrder << " interpolation on integrated histogram." << endl
3470 << " To override this choice of technique use argument ScanNone(), to change scan parameters use ScanParameters(nbins,order) argument" << endl ;
3471 }
3472
3473 return isNum ? createScanCdf(iset,nset,numScanBins,intOrder) : createIntRI(iset,nset) ;
3474 }
3475 return 0 ;
3476}
3477
3478RooAbsReal* RooAbsPdf::createScanCdf(const RooArgSet& iset, const RooArgSet& nset, Int_t numScanBins, Int_t intOrder)
3479{
3480 string name = string(GetName()) + "_NUMCDF_" + integralNameSuffix(iset,&nset).Data() ;
3481 RooRealVar* ivar = (RooRealVar*) iset.first() ;
3482 ivar->setBins(numScanBins,"numcdf") ;
3483 RooNumCdf* ret = new RooNumCdf(name.c_str(),name.c_str(),*this,*ivar,"numcdf") ;
3484 ret->setInterpolationOrder(intOrder) ;
3485 return ret ;
3486}
3487
3488
3489
3490
3491////////////////////////////////////////////////////////////////////////////////
3492/// This helper function finds and collects all constraints terms of all component p.d.f.s
3493/// and returns a RooArgSet with all those terms.
3494
3495RooArgSet* RooAbsPdf::getAllConstraints(const RooArgSet& observables, RooArgSet& constrainedParams, Bool_t stripDisconnected) const
3496{
3497 RooArgSet* ret = new RooArgSet("AllConstraints") ;
3498
3499 std::unique_ptr<RooArgSet> comps(getComponents());
3500 for (const auto arg : *comps) {
3501 auto pdf = dynamic_cast<const RooAbsPdf*>(arg) ;
3502 if (pdf && !ret->find(pdf->GetName())) {
3503 std::unique_ptr<RooArgSet> compRet(pdf->getConstraints(observables,constrainedParams,stripDisconnected));
3504 if (compRet) {
3505 ret->add(*compRet,kFALSE) ;
3506 }
3507 }
3508 }
3509
3510 return ret ;
3511}
3512
3513
3514////////////////////////////////////////////////////////////////////////////////
3515/// Returns the default numeric MC generator configuration for all RooAbsReals
3516
3518{
3520}
3521
3522
3523////////////////////////////////////////////////////////////////////////////////
3524/// Returns the specialized integrator configuration for _this_ RooAbsReal.
3525/// If this object has no specialized configuration, a null pointer is returned
3526
3528{
3529 return _specGeneratorConfig ;
3530}
3531
3532
3533
3534////////////////////////////////////////////////////////////////////////////////
3535/// Returns the specialized integrator configuration for _this_ RooAbsReal.
3536/// If this object has no specialized configuration, a null pointer is returned,
3537/// unless createOnTheFly is kTRUE in which case a clone of the default integrator
3538/// configuration is created, installed as specialized configuration, and returned
3539
3541{
3542 if (!_specGeneratorConfig && createOnTheFly) {
3544 }
3545 return _specGeneratorConfig ;
3546}
3547
3548
3549
3550////////////////////////////////////////////////////////////////////////////////
3551/// Return the numeric MC generator configuration used for this object. If
3552/// a specialized configuration was associated with this object, that configuration
3553/// is returned, otherwise the default configuration for all RooAbsReals is returned
3554
3556{
3557 const RooNumGenConfig* config = specialGeneratorConfig() ;
3558 if (config) return config ;
3559 return defaultGeneratorConfig() ;
3560}
3561
3562
3563
3564////////////////////////////////////////////////////////////////////////////////
3565/// Set the given configuration as default numeric MC generator
3566/// configuration for this object
3567
3569{
3571 delete _specGeneratorConfig ;
3572 }
3573 _specGeneratorConfig = new RooNumGenConfig(config) ;
3574}
3575
3576
3577
3578////////////////////////////////////////////////////////////////////////////////
3579/// Remove the specialized numeric MC generator configuration associated
3580/// with this object
3581
3583{
3585 delete _specGeneratorConfig ;
3586 }
3588}
3589
3590
3591
3592////////////////////////////////////////////////////////////////////////////////
3593
3595{
3596 delete _genContext ;
3597}
3598
3599
3600////////////////////////////////////////////////////////////////////////////////
3601
3602RooAbsPdf::GenSpec::GenSpec(RooAbsGenContext* context, const RooArgSet& whatVars, RooDataSet* protoData, Int_t nGen,
3603 Bool_t extended, Bool_t randProto, Bool_t resampleProto, TString dsetName, Bool_t init) :
3604 _genContext(context), _whatVars(whatVars), _protoData(protoData), _nGen(nGen), _extended(extended),
3605 _randProto(randProto), _resampleProto(resampleProto), _dsetName(dsetName), _init(init)
3606{
3607}
3608
3609
3610
3611////////////////////////////////////////////////////////////////////////////////
3612
3613void RooAbsPdf::setNormRange(const char* rangeName)
3614{
3615 if (rangeName) {
3616 _normRange = rangeName ;
3617 } else {
3618 _normRange.Clear() ;
3619 }
3620
3621 if (_norm) {
3623 _norm = 0 ;
3624 }
3625}
3626
3627
3628////////////////////////////////////////////////////////////////////////////////
3629
3630void RooAbsPdf::setNormRangeOverride(const char* rangeName)
3631{
3632 if (rangeName) {
3633 _normRangeOverride = rangeName ;
3634 } else {
3636 }
3637
3638 if (_norm) {
3640 _norm = 0 ;
3641 }
3642}
header file containing the templated implementation of matrix inversion routines for use with ROOT's ...
void Class()
Definition: Class.C:29
#define e(i)
Definition: RSha256.hxx:103
#define coutI(a)
Definition: RooMsgService.h:31
#define cxcoutI(a)
Definition: RooMsgService.h:86
#define cxcoutD(a)
Definition: RooMsgService.h:82
#define oocoutW(o, a)
Definition: RooMsgService.h:48
#define coutW(a)
Definition: RooMsgService.h:33
#define coutE(a)
Definition: RooMsgService.h:34
#define ccoutI(a)
Definition: RooMsgService.h:39
#define ccoutD(a)
Definition: RooMsgService.h:38
double _rf_fast_log(double x)
Definition: RooVDTHeaders.h:51
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
double sqrt(double)
double log(double)
char * Form(const char *fmt,...)
Status_t status(std::size_t begin, std::size_t size) const
Return the status of the batch starting at begin.
Definition: BatchData.h:48
RooSpan< const double > getBatch(std::size_t begin, std::size_t batchSize) const
Definition: BatchData.cxx:25
bool setStatus(std::size_t begin, std::size_t size, Status_t stat)
Set the status of a batch with the given start point and size.
Definition: BatchData.h:84
void markDirty()
Mark all batches dirty. This will trigger recomputations.
Definition: BatchData.h:97
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
class to compute the Cholesky decomposition of a matrix
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:240
void clearValueAndShapeDirty() const
Definition: RooAbsArg.h:521
friend class RooDataSet
Definition: RooAbsArg.h:599
RooWorkspace * _myws
Prevent 'AlwaysDirty' mode for this node.
Definition: RooAbsArg.h:653
Bool_t dependsOn(const RooAbsCollection &serverList, const RooAbsArg *ignoreArg=0, Bool_t valueOnly=kFALSE) const
Test whether we depend on (ie, are served by) any object in the specified collection.
Definition: RooAbsArg.cxx:734
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:83
friend class RooArgSet
Definition: RooAbsArg.h:551
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:281
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
Definition: RooAbsArg.cxx:306
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:548
Bool_t isValueDirty() const
Definition: RooAbsArg.h:390
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1917
virtual void constOptimizeTestStatistic(ConstOpCode opcode, Bool_t doAlsoTrackingOpt=kTRUE)
Interface function signaling a request to perform constant term optimization.
Definition: RooAbsArg.cxx:1713
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Change cache operation mode to given mode.
Definition: RooAbsArg.cxx:1726
OperMode _operMode
Mark batches as dirty (only meaningful for RooAbsReal).
Definition: RooAbsArg.h:635
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2115
bool _allBatchesDirty
Definition: RooAbsArg.h:633
Bool_t getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
Definition: RooAbsArg.cxx:284
void setProxyNormSet(const RooArgSet *nset)
Forward a change in the cached normalization argset to all the registered proxies.
Definition: RooAbsArg.cxx:1260
friend class RooProjectedPdf
Definition: RooAbsArg.h:558
RefCountList_t _serverList
Definition: RooAbsArg.h:559
RooArgSet * getComponents() const
Definition: RooAbsArg.cxx:678
Bool_t isConstant() const
Definition: RooAbsArg.h:320
void branchNodeServerList(RooAbsCollection *list, const RooAbsArg *arg=0, Bool_t recurseNonDerived=kFALSE) const
Fill supplied list with all branch nodes of the arg tree starting with ourself as top node.
Definition: RooAbsArg.cxx:490
RooAbsArg * findServer(const char *name) const
Return server of this with name name. Returns nullptr if not found.
Definition: RooAbsArg.h:167
OperMode operMode() const
Definition: RooAbsArg.h:453
void setInterpolationOrder(Int_t order)
Set interpolation order of RooHistFunct representing cache histogram.
RooAbsCategoryLValue is the common abstract base class for objects that represent a discrete value th...
RooAbsCollection * selectCommon(const RooAbsCollection &refColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
Int_t getSize() const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
RooAbsArg * first() const
void reserve(Storage_t::size_type count)
RooAbsCollection * selectByName(const char *nameList, Bool_t verbose=kFALSE) const
Create a subset of the current collection, consisting only of those elements with names matching the ...
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
RooAbsCollection * selectByAttrib(const char *name, Bool_t value) const
Create a subset of the current collection, consisting only of those elements with the specified attri...
std::string contentsString() const
Return comma separated list of contained object names as STL string.
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:39
virtual const RooArgSet * get() const
Definition: RooAbsData.h:82
virtual Bool_t isNonPoissonWeighted() const
Definition: RooAbsData.h:103
virtual Double_t weight() const =0
virtual Bool_t isWeighted() const
Definition: RooAbsData.h:99
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:307
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
Bool_t isValid() const
virtual RooDataSet * generate(Double_t nEvents=0, Bool_t skipInit=kFALSE, Bool_t extendedMode=kFALSE)
Generate the specified number of events with nEvents>0 and and return a dataset containing the genera...
virtual void setExpectedData(Bool_t)
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of prototype data to that in the lookup tables passed as argument.
Normalization set with for above integral.
Definition: RooAbsPdf.h:324
virtual ~CacheElem()
Destructor of normalization cache element.
Definition: RooAbsPdf.cxx:3344
void operModeHook(RooAbsArg::OperMode)
Dummy implementation.
Definition: RooAbsPdf.cxx:3333
RooAbsReal * _norm
Definition: RooAbsPdf.h:330
RooArgSet _whatVars
Definition: RooAbsPdf.h:81
RooAbsGenContext * _genContext
Definition: RooAbsPdf.h:80
Bool_t _resampleProto
Definition: RooAbsPdf.h:86
RooDataSet * _protoData
Definition: RooAbsPdf.h:82
GenSpec * prepareMultiGen(const RooArgSet &whatVars, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none())
Prepare GenSpec configuration object for efficient generation of multiple datasets from identical spe...
Definition: RooAbsPdf.cxx:2195
RooSpan< const double > getLogValBatch(std::size_t begin, std::size_t batchSize, const RooArgSet *normSet=nullptr) const
Compute the log-likelihoods for all events in the requested batch.
Definition: RooAbsPdf.cxx:781
RooObjCacheManager _normMgr
Definition: RooAbsPdf.h:332
Double_t getNorm(const RooArgSet &nset) const
Definition: RooAbsPdf.h:208
friend class CacheElem
Definition: RooAbsPdf.h:334
virtual void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
Definition: RooAbsPdf.cxx:2453
virtual ~RooAbsPdf()
Destructor.
Definition: RooAbsPdf.cxx:263
void logBatchComputationErrors(RooSpan< const double > &outputs, std::size_t begin) const
Scan through outputs and fix+log all nans and negative values.
Definition: RooAbsPdf.cxx:758
void setGeneratorConfig()
Remove the specialized numeric MC generator configuration associated with this object.
Definition: RooAbsPdf.cxx:3582
virtual void resetErrorCounters(Int_t resetValue=10)
Reset error counter to given value, limiting the number of future error messages for this pdf to 'res...
Definition: RooAbsPdf.cxx:662
virtual RooAbsReal * createChi2(RooDataHist &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Create a from a histogram and this function.
Definition: RooAbsPdf.cxx:1882
virtual RooFitResult * chi2FitTo(RooDataHist &data, const RooLinkedList &cmdList)
Calls RooAbsPdf::createChi2(RooDataSet& data, const RooLinkedList& cmdList) and returns fit result.
Definition: RooAbsPdf.cxx:1837
static int verboseEval()
Return global level of verbosity for p.d.f. evaluations.
Definition: RooAbsPdf.cxx:3323
Bool_t traceEvalPdf(Double_t value) const
Check that passed value is positive and not 'not-a-number'.
Definition: RooAbsPdf.cxx:437
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:957
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const
Interface function to create a generator context from a p.d.f.
Definition: RooAbsPdf.cxx:2033
RooAbsReal * createScanCdf(const RooArgSet &iset, const RooArgSet &nset, Int_t numScanBins, Int_t intOrder)
Definition: RooAbsPdf.cxx:3478
TString _normRange
MC generator configuration specific for this object.
Definition: RooAbsPdf.h:355
void setNormRange(const char *rangeName)
Definition: RooAbsPdf.cxx:3613
virtual RooDataHist * generateBinned(const RooArgSet &whatVars, Double_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none()) const
As RooAbsPdf::generateBinned(const RooArgSet&, const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,...
Definition: RooAbsPdf.h:104
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooAbsPdf.cxx:2431
RooNumGenConfig * specialGeneratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
Definition: RooAbsPdf.cxx:3527
virtual RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize, const RooArgSet *normSet=nullptr) const
Compute batch of values for given range, and normalise by integrating over the observables in nset.
Definition: RooAbsPdf.cxx:348
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:237
virtual Bool_t selfNormalized() const
Definition: RooAbsPdf.h:223
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Print multi line detailed information of this RooAbsPdf.
Definition: RooAbsPdf.cxx:2006
Double_t analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const
Analytical integral with normalization (see RooAbsReal::analyticalIntegralWN() for further informatio...
Definition: RooAbsPdf.cxx:417
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition: RooAbsPdf.h:55
virtual RooAbsGenContext * autoGenContext(const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE, Bool_t autoBinned=kTRUE, const char *binnedTag="") const
Definition: RooAbsPdf.cxx:2042
Int_t _traceCount
Definition: RooAbsPdf.h:348
RooAbsReal * _norm
Definition: RooAbsPdf.h:321
virtual Double_t extendedTerm(Double_t observedEvents, const RooArgSet *nset=0) const
Return the extended likelihood term ( ) of this PDF for the given number of observed events.
Definition: RooAbsPdf.cxx:817
virtual void printValue(std::ostream &os) const
Print value of p.d.f, also print normalization integral that was last used, if any.
Definition: RooAbsPdf.cxx:1990
virtual RooArgSet * getAllConstraints(const RooArgSet &observables, RooArgSet &constrainedParams, Bool_t stripDisconnected=kTRUE) const
This helper function finds and collects all constraints terms of all component p.d....
Definition: RooAbsPdf.cxx:3495
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1286
virtual Bool_t syncNormalization(const RooArgSet *dset, Bool_t adjustProxies=kTRUE) const
Verify that the normalization integral cached with this PDF is valid for given set of normalization o...
Definition: RooAbsPdf.cxx:530
virtual Double_t getValV(const RooArgSet *set=0) const
Return current value, normalized by integrating over the observables in nset.
Definition: RooAbsPdf.cxx:281
Int_t _errorCount
Definition: RooAbsPdf.h:344
Bool_t _selectComp
Definition: RooAbsPdf.h:351
@ CanBeExtended
Definition: RooAbsPdf.h:230
@ MustBeExtended
Definition: RooAbsPdf.h:230
@ CanNotBeExtended
Definition: RooAbsPdf.h:230
virtual RooPlot * paramOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Add a box with parameter values (and errors) to the specified frame.
Definition: RooAbsPdf.cxx:3152
void setTraceCounter(Int_t value, Bool_t allNodes=kFALSE)
Reset trace counter to given value, limiting the number of future trace messages for this pdf to 'val...
Definition: RooAbsPdf.cxx:674
RooAbsReal * createCdf(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Create a cumulative distribution function of this p.d.f in terms of the observables listed in iset.
Definition: RooAbsPdf.cxx:3400
virtual const RooAbsReal * getNormObj(const RooArgSet *set, const RooArgSet *iset, const TNamed *rangeName=0) const
Return pointer to RooAbsReal object that implements calculation of integral over observables iset in ...
Definition: RooAbsPdf.cxx:495
virtual Double_t getLogVal(const RooArgSet *set=0) const
Return the log of the current value with given normalization An error message is printed if the argum...
Definition: RooAbsPdf.cxx:700
Int_t * randomizeProtoOrder(Int_t nProto, Int_t nGen, Bool_t resample=kFALSE) const
Return lookup table with randomized access order for prototype events, given nProto prototype data ev...
Definition: RooAbsPdf.cxx:2387
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:3303
Int_t _negCount
Definition: RooAbsPdf.h:349
void setNormRangeOverride(const char *rangeName)
Definition: RooAbsPdf.cxx:3630
virtual RooDataSet * generateSimGlobal(const RooArgSet &whatVars, Int_t nEvents)
Special generator interface for generation of 'global observables' – for RooStats tools.
Definition: RooAbsPdf.cxx:2689
virtual Bool_t isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
Definition: RooAbsPdf.cxx:2466
const RooNumGenConfig * getGeneratorConfig() const
Return the numeric MC generator configuration used for this object.
Definition: RooAbsPdf.cxx:3555
virtual void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
Definition: RooAbsPdf.cxx:2441
virtual ExtendMode extendMode() const
Definition: RooAbsPdf.h:231
RooAbsPdf()
Default constructor.
Definition: RooAbsPdf.cxx:205
virtual RooAbsGenContext * binnedGenContext(const RooArgSet &vars, Bool_t verbose=kFALSE) const
Return a binned generator context.
Definition: RooAbsPdf.cxx:2023
static RooNumGenConfig * defaultGeneratorConfig()
Returns the default numeric MC generator configuration for all RooAbsReals.
Definition: RooAbsPdf.cxx:3517
RooNumGenConfig * _specGeneratorConfig
Definition: RooAbsPdf.h:353
RooArgSet * _normSet
Normalization integral (owned by _normMgr)
Definition: RooAbsPdf.h:322
Double_t _rawValue
Definition: RooAbsPdf.h:320
static TString _normRangeOverride
Definition: RooAbsPdf.h:356
static Int_t _verboseEval
Definition: RooAbsPdf.h:315
virtual RooAbsPdf * createProjection(const RooArgSet &iset)
Return a p.d.f that represent a projection of this p.d.f integrated over given observables.
Definition: RooAbsPdf.cxx:3362
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition: RooAbsPdf.h:118
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Bool_t hasRange(const char *name) const
Check if variable has a binning with given name.
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
void plotOnCompSelect(RooArgSet *selNodes) const
Helper function for plotting of composite p.d.fs.
@ RelativeExpected
Definition: RooAbsReal.h:249
TString getTitle(Bool_t appendUnit=kFALSE) const
Return this variable's title string.
Definition: RooAbsReal.cxx:254
RooDerivative * derivative(RooRealVar &obs, Int_t order=1, Double_t eps=0.001)
Return function representing first, second or third order derivative of this function.
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg(), const RooCmdArg &arg10=RooCmdArg()) const
Plot (project) PDF on specified frame.
virtual Double_t evaluate() const =0
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
TString integralNameSuffix(const RooArgSet &iset, const RooArgSet *nset=0, const char *rangeName=0, Bool_t omitEmpty=kFALSE) const
Construct string with unique suffix name to give to integral object that encodes integrated observabl...
Definition: RooAbsReal.cxx:809
virtual RooSpan< double > evaluateBatch(std::size_t begin, std::size_t maxSize) const
Evaluate function for a batch of input data points.
RooAbsReal * createIntRI(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Utility function for createRunningIntegral.
RooFitResult * chi2FitDriver(RooAbsReal &fcn, RooLinkedList &cmdList)
Internal driver function for chi2 fits.
Double_t _value
Definition: RooAbsReal.h:446
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Structure printing.
Definition: RooAbsReal.cxx:477
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
virtual void enableOffsetting(Bool_t)
Definition: RooAbsReal.h:339
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:562
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
Bool_t plotSanityChecks(RooPlot *frame) const
Utility function for plotOn(), perform general sanity check on frame to ensure safe plotting operatio...
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooAbsReal.cxx:416
const RooNumIntConfig * getIntegratorConfig() const
Return the numeric integration configuration used for this object.
virtual Bool_t isBinnedDistribution(const RooArgSet &) const
Definition: RooAbsReal.h:309
void logEvalError(const char *message, const char *serverValueString=0) const
Log evaluation error message.
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, Double_t scaleFactor, Bool_t correctForBinVolume=kFALSE, Bool_t showProgress=kFALSE) const
Fill a RooDataHist with values sampled from this function at the bin centers.
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
Definition: RooAddition.h:26
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooBinnedGenContext is an efficient implementation of the generator context specific for binned pdfs.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
RooCachedReal is an implementation of RooAbsCachedReal that can cache any external RooAbsReal input f...
Definition: RooCachedReal.h:20
void setCacheSource(Bool_t flag)
Definition: RooCachedReal.h:44
Class RooChi2Var implements a simple chi^2 calculation from a binned dataset and a PDF.
Definition: RooChi2Var.h:25
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:28
void setInt(Int_t idx, Int_t value)
Definition: RooCmdArg.h:67
void setString(Int_t idx, const char *value)
Definition: RooCmdArg.h:73
Class RooCmdConfig is a configurable parser for RooCmdArg named arguments.
Definition: RooCmdConfig.h:27
RooConstraintSum calculates the sum of the -(log) likelihoods of a set of RooAbsPfs that represent co...
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual Double_t weight() const
Definition: RooDataHist.h:106
void set(Double_t weight, Double_t wgtErr=-1)
Set the weight and weight error of the bin enclosing the current (i.e.
virtual Double_t sumEntries() const
virtual Int_t numEntries() const
Return the number of bins.
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
void SetName(const char *name)
Change the name of the RooDataHist.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
void SetName(const char *name) override
Change the name of this dataset into the given name.
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
Definition: RooGenContext.h:30
RooInt is a minimal implementation of a TObject holding a Int_t value.
Definition: RooInt.h:22
RooLinkedList is an collection class for internal use, storing a collection of RooAbsArg pointers in ...
Definition: RooLinkedList.h:36
TObject * FindObject(const char *name) const
Return pointer to obejct with given name.
virtual void Add(TObject *arg)
Definition: RooLinkedList.h:63
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:38
RooMinuit is a wrapper class around TFitter/TMinuit that provides a seamless interface between the MI...
Definition: RooMinuit.h:39
Class RooNLLVar implements a -log(likelihood) calculation from a dataset and a PDF.
Definition: RooNLLVar.h:26
static const char * str(const TNamed *ptr)
Return C++ string corresponding to given TNamed pointer.
Definition: RooNameReg.cxx:103
RooNameSet is a utility class that stores the names the objects in a RooArget.
Definition: RooNameSet.h:24
RooArgSet * select(const RooArgSet &list) const
Construct a RooArgSet of objects in input 'list' whose names match to those in the internal name list...
Definition: RooNameSet.cxx:187
void setNameList(const char *givenList)
Definition: RooNameSet.cxx:144
const char * content() const
Definition: RooNameSet.h:50
Class RooNumCdf is an implementation of RooNumRunningInt specialized to calculate cumulative distribu...
Definition: RooNumCdf.h:17
RooNumGenConfig holds the configuration parameters of the various numeric integrators used by RooReal...
static RooNumGenConfig & defaultConfig()
Return reference to instance of default numeric integrator configuration object.
void sterilize()
Clear the cache payload but retain slot mapping w.r.t to normalization and integration sets.
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
void addObject(TObject *obj, Option_t *drawOptions="", Bool_t invisible=kFALSE)
Add a generic object to this plot.
Definition: RooPlot.cxx:443
Double_t getFitRangeNEvt() const
Return the number of events in the fit range.
Definition: RooPlot.h:141
Double_t getFitRangeBinW() const
Return the bin width that is being used to normalise the PDF.
Definition: RooPlot.h:144
const RooArgSet * getNormVars() const
Definition: RooPlot.h:148
RooAbsRealLValue * getPlotVar() const
Definition: RooPlot.h:139
void updateNormVars(const RooArgSet &vars)
Install the given set of observables are reference normalization variables for this frame.
Definition: RooPlot.cxx:401
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
static UInt_t integer(UInt_t max, TRandom *generator=randomGenerator())
Return an integer uniformly distributed from [0,n-1].
Definition: RooRandom.cxx:102
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
RooRealIntegral performs hybrid numerical/analytical integrals of RooAbsReal objects.
const RooArgSet & numIntRealVars() const
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:430
void setRange(const char *name, Double_t min, Double_t max)
Set range named 'name to [min,max].
Definition: RooRealVar.cxx:539
TString * format(const RooCmdArg &formatArg) const
Format contents of RooRealVar for pretty printing on RooPlot parameter boxes.
Definition: RooRealVar.cxx:870
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
constexpr std::span< T >::index_type size() const noexcept
Definition: RooSpan.h:123
Bool_t defineSetInternal(const char *name, const RooArgSet &aset)
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
TString fName
Definition: TNamed.h:32
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
A Pave (see TPave) with text, lines or/and boxes inside.
Definition: TPaveText.h:21
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
virtual UInt_t Integer(UInt_t imax)
Returns a random integer uniformly distributed on the interval [ 0, imax-1 ].
Definition: TRandom.cxx:349
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
void Clear()
Clear string without changing its capacity.
Definition: TString.cxx:1176
const char * Data() const
Definition: TString.h:364
TString & Append(const char *cs)
Definition: TString.h:559
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
EvaluateInfo init(std::vector< RooRealProxy > parameters, std::vector< ArrayWrapper * > wrappers, std::vector< double * > arrays, size_t begin, size_t batchSize)
double T(double x)
Definition: ChebyshevPol.h:34
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
RooCmdArg WeightVar(const char *name, Bool_t reinterpretAsWeight=kFALSE)
RooCmdArg NormRange(const char *rangeNameList)
RooCmdArg Hesse(Bool_t flag=kTRUE)
@ Minimization
Definition: RooGlobalFunc.h:67
@ Generation
Definition: RooGlobalFunc.h:67
@ NumIntegration
Definition: RooGlobalFunc.h:69
@ InputArguments
Definition: RooGlobalFunc.h:68
RooCmdArg SupNormSet(const RooArgSet &nset)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
RooCmdArg Normalization(Double_t scaleFactor)
std::vector< std::string > tokenise(const std::string &str, const std::string &delims)
Tokenise the string by splitting at the characters in delims.
Definition: RooHelpers.cxx:25
static constexpr double pc
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
Definition: first.py:1
Double_t scaleFactor
Definition: RooAbsReal.h:475
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
static void output(int code)
Definition: gifencode.c:226