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