Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooSimultaneous.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\file RooSimultaneous.cxx
19\class RooSimultaneous
20\ingroup Roofitcore
21
22Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
23The class takes an index category, which is used as a selector
24for PDFs, and a list of PDFs, each associated
25with a state of the index category. RooSimultaneous always returns
26the value of the PDF that is associated with the current value
27of the index category.
28
29Extended likelihood fitting is supported if all components support
30extended likelihood mode. If the index category is among the normalization
31variables, the expected number of events by a RooSimultaneous is the sum over
32all components. Otherwise, it is that of the component p.d.f. selected by the
33index category. For a RooSuperCategory index, it counts as being among the
34normalization variables if all of its input categories are.
35
36The index category can be accessed using indexCategory().
37
38###Generating events
39When generating events from a RooSimultaneous, the index category has to be added to
40the dataset. Further, the PDF needs to know the relative probabilities of each category, i.e.,
41how many events are in which category. This can be achieved in two ways:
42- Generating with proto data that have category entries: An event from the same category as
43in the proto data is created for each event in the proto data.
44See RooAbsPdf::generate(const RooArgSet&,const RooDataSet&,Int_t,bool,bool,bool) const.
45- No proto data: A category is chosen randomly.
46\note This requires that the PDFs building the simultaneous are extended. In this way,
47the relative probability of each category can be calculated from the number of events
48in each category.
49
50###Switch mode
51If the index category is not an observable, but a parameter, the RooSimultaneous
52acts as a "switch" that evaluates to the component selected by the current index
53state, like a RooMultiPdf. Its expected number of events is then the one of the
54selected component.
55
56###Nesting simultaneous PDFs
57RooSimultaneous components are only supported if their index category is a
58parameter (switch mode); otherwise, creating the likelihood fails. To combine
59channels, use one RooSimultaneous over a RooCategory with one state per channel,
60or over a RooSuperCategory if the channels span the full product of several
61categories. RooSimultaneous objects flattened by ROOT 6.42 or earlier and read
62from files still work.
63**/
64
65#include "RooSimultaneous.h"
66
67#include "Roo1DTable.h"
69#include "RooAbsData.h"
70#include "RooAddPdf.h"
71#include "RooArgSet.h"
72#include "RooBinSamplingPdf.h"
73#include "RooCategory.h"
74#include "RooCmdConfig.h"
76#include "RooDataHist.h"
77#include "RooDataSet.h"
78#include "RooFormulaVar.h"
79#include "RooGlobalFunc.h"
80#include "RooMsgService.h"
81#include "RooNameReg.h"
82#include "RooPlot.h"
83#include "RooRandom.h"
84#include "RooRealVar.h"
85#include "RooSimGenContext.h"
87#include "RooSuperCategory.h"
88
89#include "RooFitImplHelpers.h"
90
91#include <ROOT/StringUtils.hxx>
92
93#include <algorithm>
94#include <iostream>
95#include <map>
96
97namespace {
98
99std::map<std::string, RooAbsPdf *> createPdfMap(const RooArgList &inPdfList, RooAbsCategoryLValue &inIndexCat)
100{
101 std::map<std::string, RooAbsPdf *> pdfMap;
102 auto indexCatIt = inIndexCat.begin();
103 for (unsigned int i = 0; i < inPdfList.size(); ++i) {
104 auto pdf = static_cast<RooAbsPdf *>(&inPdfList[i]);
105 const auto &nameIdx = (*indexCatIt++);
106 pdfMap[nameIdx.first] = pdf;
107 }
108 return pdfMap;
109}
110
111void replaceOrAdd(RooLinkedList &lst, TObject &obj)
112{
113 TObject *old = lst.FindObject(obj.GetName());
114 if (old)
115 lst.Replace(old, &obj);
116 else
117 lst.Add(&obj);
118}
119
120} // namespace
121
122using std::string;
123
124
125
126
127////////////////////////////////////////////////////////////////////////////////
128/// Constructor with index category. PDFs associated with indexCat
129/// states can be added after construction with the addPdf() function.
130///
131/// RooSimultaneous can function without having a PDF associated
132/// with every single state. The normalization in such cases is taken
133/// from the number of registered PDFs, but getVal() will assert if
134/// when called for an unregistered index state.
135
136RooSimultaneous::RooSimultaneous(const char *name, const char *title,
137 RooAbsCategoryLValue& inIndexCat) :
138 RooSimultaneous{name, title, std::map<std::string, RooAbsPdf*>{}, inIndexCat}
139{
140}
141
142
143////////////////////////////////////////////////////////////////////////////////
144/// Constructor from index category and full list of PDFs.
145/// In this constructor form, a PDF must be supplied for each indexCat state
146/// to avoid ambiguities. The PDFs are associated with the states of the
147/// index category as they appear when iterating through the category states
148/// with RooAbsCategory::begin() and RooAbsCategory::end(). This usually means
149/// they are associated by ascending index numbers.
150///
151/// PDFs may not overlap (i.e. share any variables) with the index category (function)
152
153RooSimultaneous::RooSimultaneous(const char *name, const char *title,
154 const RooArgList& inPdfList, RooAbsCategoryLValue& inIndexCat) :
155 RooSimultaneous{name, title, createPdfMap(inPdfList, inIndexCat), inIndexCat}
156{
157 if (inPdfList.size() != inIndexCat.size()) {
158 std::stringstream errMsg;
159 errMsg << "RooSimultaneous::ctor(" << GetName()
160 << " ERROR: Number PDF list entries must match number of index category states, no PDFs added";
161 coutE(InputArguments) << errMsg.str() << std::endl;
162 throw std::invalid_argument(errMsg.str());
163 }
164}
165
166
167////////////////////////////////////////////////////////////////////////////////
168
169RooSimultaneous::RooSimultaneous(const char *name, const char *title, std::map<string, RooAbsPdf *> pdfMap,
170 RooAbsCategoryLValue &inIndexCat)
171 : RooAbsPdf(name, title),
172 _plotCoefNormSet("!plotCoefNormSet", "plotCoefNormSet", this, false, false),
173 _partIntMgr(this, 10),
174 _indexCat("indexCat", "Index category", this, inIndexCat)
175{
176 for (auto const &item : pdfMap) {
177 addPdf(*item.second, item.first.c_str());
178 }
179}
180
181/// For internal use in RooFit.
182RooSimultaneous::RooSimultaneous(const char *name, const char *title,
184 RooAbsCategoryLValue &inIndexCat)
185 : RooSimultaneous(name, title, RooFit::Detail::flatMapToStdMap(pdfMap), inIndexCat)
186{
187}
188
189////////////////////////////////////////////////////////////////////////////////
190/// Copy constructor
191
193 RooAbsPdf(other,name),
194 _plotCoefNormSet("!plotCoefNormSet",this,other._plotCoefNormSet),
195 _plotCoefNormRange(other._plotCoefNormRange),
196 _partIntMgr(other._partIntMgr,this),
197 _indexCat("indexCat",this,other._indexCat),
198 _numPdf(other._numPdf)
199{
200 // Copy proxy list
201 for(auto* proxy : static_range_cast<RooRealProxy*>(other._pdfProxyList)) {
202 _pdfProxyList.Add(new RooRealProxy(proxy->GetName(),this,*proxy)) ;
203 }
204}
205
206
207
208////////////////////////////////////////////////////////////////////////////////
209/// Destructor
210
215
216
217
218////////////////////////////////////////////////////////////////////////////////
219/// Return the p.d.f associated with the given index category name
220
222{
223 RooRealProxy* proxy = static_cast<RooRealProxy*>(_pdfProxyList.FindObject(catName.c_str()));
224 return proxy ? static_cast<RooAbsPdf*>(proxy->absArg()) : nullptr;
225}
226
227
228
229////////////////////////////////////////////////////////////////////////////////
230/// Associate given PDF with index category state label 'catLabel'.
231/// The name state must be already defined in the index category.
232///
233/// RooSimultaneous can function without having a PDF associated
234/// with every single state. The normalization in such cases is taken
235/// from the number of registered PDFs, but getVal() will fail if
236/// called for an unregistered index state.
237///
238/// PDFs may not overlap (i.e. share any variables) with the index category (function).
239/// \param[in] pdf PDF to be added.
240/// \param[in] catLabel Name of the category state to be associated to the PDF.
241/// \return `true` in case of failure.
242
243bool RooSimultaneous::addPdf(const RooAbsPdf& pdf, const char* catLabel)
244{
245 // PDFs cannot overlap with the index category
246 if (pdf.dependsOn(_indexCat.arg())) {
247 coutE(InputArguments) << "RooSimultaneous::addPdf(" << GetName() << "): PDF '" << pdf.GetName()
248 << "' overlaps with index category '" << _indexCat.arg().GetName() << "'."<< std::endl ;
249 return true ;
250 }
251
252 // Each index state can only have one PDF associated with it
253 if (_pdfProxyList.FindObject(catLabel)) {
254 coutE(InputArguments) << "RooSimultaneous::addPdf(" << GetName() << "): index state '"
255 << catLabel << "' has already an associated PDF." << std::endl ;
256 return true ;
257 }
258
259 // Create a proxy named after the associated index state
260 TObject* proxy = new RooRealProxy(catLabel,catLabel,this,const_cast<RooAbsPdf&>(pdf));
261 _pdfProxyList.Add(proxy) ;
262 _numPdf += 1 ;
263
264 return false ;
265}
266
267////////////////////////////////////////////////////////////////////////////////
268/// Examine the pdf components and check if one of them can be extended or must be extended.
269/// It is enough to have one component that can be extended or must be extended to return the flag in
270/// the total simultaneous pdf.
271
273{
274 bool anyCanExtend = false;
275
276 for (auto *proxy : static_range_cast<RooRealProxy *>(_pdfProxyList)) {
277 auto &pdf = static_cast<RooAbsPdf const&>(proxy->arg());
278 if (pdf.mustBeExtended())
279 return MustBeExtended;
280 anyCanExtend |= pdf.canBeExtended();
281 }
282 return anyCanExtend ? CanBeExtended : CanNotBeExtended;
283}
284
285////////////////////////////////////////////////////////////////////////////////
286/// Return the current value:
287/// the value of the PDF associated with the current index category state
288
290{
291 // Retrieve the proxy by index name
293 if(!proxy) {
294 return 0;
295 }
296
297 double nEvtTot = 1.0;
298 double nEvtCat = 1.0;
299
300 // Relative yield weight only if normalized over the index states. In
301 // "switch mode", the value is the one of the selected component.
303
304 nEvtTot = 0;
305 nEvtCat = 0;
306
307 for (auto *proxy2 : static_range_cast<RooRealProxy *>(_pdfProxyList)) {
308 auto &pdf2 = static_cast<RooAbsPdf const &>(proxy2->arg());
309 if(!pdf2.canBeExtended()) {
310 // If one of the pdfs can't be expected, reset the normalization
311 // factor to one and break out of the loop.
312 nEvtTot = 1.0;
313 nEvtCat = 1.0;
314 break;
315 }
316 const double nEvt = pdf2.expectedEvents(_normSet);
317 nEvtTot += nEvt;
318 if (proxy == proxy2) {
319 // Matching by proxy by pointer rather than pdfs, because it's
320 // possible to have the same pdf used in different states.
321 nEvtCat += nEvt;
322 }
323 }
324 }
325 double catFrac = nEvtCat / nEvtTot;
326
327 // Return the selected PDF value, normalized by the relative number of
328 // expected events if applicable.
329 return *proxy * catFrac;
330}
331
332////////////////////////////////////////////////////////////////////////////////
333/// Return the number of expected events: If the index is in nset,
334/// then return the sum of the expected events of all components,
335/// otherwise return the number of expected events of the PDF
336/// associated with the current index category state
337
339{
340 if (indexCatIsInNormSet(nset)) {
341
342 double sum(0) ;
343
344 for(auto * proxy : static_range_cast<RooRealProxy*>(_pdfProxyList)) {
345 sum += (static_cast<RooAbsPdf*>(proxy->absArg()))->expectedEvents(nset) ;
346 }
347
348 return sum ;
349
350 } else {
351
352 // Retrieve the proxy by index name
354
355 //assert(proxy!=0) ;
356 if (proxy==nullptr) return 0 ;
357
358 // Return the selected PDF value, normalized by the number of index states
359 return (static_cast<RooAbsPdf*>(proxy->absArg()))->expectedEvents(nset);
360 }
361}
362
363
364
365////////////////////////////////////////////////////////////////////////////////
366/// Expected events function for the new evaluation backends: the sum over all
367/// components, or the selected component in "switch mode" (see expectedEvents()).
368
369std::unique_ptr<RooAbsReal> RooSimultaneous::createExpectedEventsFunc(const RooArgSet *nset) const
370{
371 const std::string name = std::string(GetName()) + "_expectedEvents";
372
373 // One function per distinct component, because the same pdf can be used in
374 // several states. Components that can't be extended contribute zero, like
375 // in expectedEvents().
376 RooArgList formulaArgs;
377 RooArgSet ownedFuncs;
378 std::map<RooAbsPdf const *, std::string> funcRefs;
379 bool failed = false;
380 auto ref = [&](RooAbsPdf const &pdf) -> std::string {
381 if (!pdf.canBeExtended()) {
382 return "0.0";
383 }
384 auto found = funcRefs.find(&pdf);
385 if (found != funcRefs.end()) {
386 return found->second;
387 }
388 std::unique_ptr<RooAbsReal> func = pdf.createExpectedEventsFunc(nset);
389 if (!func) {
390 failed = true;
391 return "0.0";
392 }
393 // Different components can have the same name
394 func->SetName((std::string(func->GetName()) + "_" + std::to_string(funcRefs.size())).c_str());
395 std::string const funcRef = "x[" + std::to_string(formulaArgs.size()) + "]";
396 formulaArgs.add(*func);
397 ownedFuncs.addOwned(std::move(func));
398 return funcRefs[&pdf] = funcRef;
399 };
400
401 std::string formula;
402 if (indexCatIsInNormSet(nset)) {
403 for (auto *proxy : static_range_cast<RooRealProxy *>(_pdfProxyList)) {
404 formula += (formula.empty() ? "" : " + ") + ref(static_cast<RooAbsPdf const &>(proxy->arg()));
405 }
406 if (formula.empty()) {
407 formula = "0.0";
408 }
409 } else {
410 // Nested ternary expression that selects the component matching the
411 // index state, keyed on the state index numbers.
412 formulaArgs.add(_indexCat.arg());
413 std::size_t numStates = 0;
414 for (auto const &nameIdx : _indexCat.arg()) {
415 auto *proxy = static_cast<RooRealProxy *>(_pdfProxyList.FindObject(nameIdx.first.c_str()));
416 if (!proxy) {
417 continue;
418 }
419 formula += "x[0] == " + std::to_string(nameIdx.second) + " ? " +
420 ref(static_cast<RooAbsPdf const &>(proxy->arg())) + " : (";
421 ++numStates;
422 }
423 formula += "0.0" + std::string(numStates, ')');
424 }
425
426 // Like the default implementation, return nullptr if a component can't
427 // create the function.
428 if (failed) {
429 return nullptr;
430 }
431
432 auto out = std::make_unique<RooFormulaVar>(name.c_str(), name.c_str(), formula.c_str(), formulaArgs);
433 out->addOwnedComponents(std::move(ownedFuncs));
434 return out;
435}
436
437////////////////////////////////////////////////////////////////////////////////
438/// Forward determination of analytical integration capabilities to component p.d.f.s
439/// A unique code is assigned to the combined integration capabilities of all associated
440/// p.d.f.s
441
443 const RooArgSet* normSet, const char* rangeName) const
444{
445 // Declare that we can analytically integrate all requested observables
446 analVars.add(allVars) ;
447
448 // Retrieve (or create) the required partial integral list
449 Int_t code ;
450
451 // Check if this configuration was created before
452 CacheElem* cache = static_cast<CacheElem*>(_partIntMgr.getObj(normSet,&analVars,nullptr,RooNameReg::ptr(rangeName))) ;
453 if (cache) {
454 code = _partIntMgr.lastIndex() ;
455 return code+1 ;
456 }
457 cache = new CacheElem ;
458
459 // Create the partial integral set for this request
460 for(auto * proxy : static_range_cast<RooRealProxy*>(_pdfProxyList)) {
461 cache->_partIntList.addOwned(std::unique_ptr<RooAbsReal>{proxy->arg().createIntegral(analVars,normSet,nullptr,rangeName)});
462 }
463
464 // Store the partial integral list and return the assigned code ;
465 code = _partIntMgr.setObj(normSet,&analVars,cache,RooNameReg::ptr(rangeName)) ;
466
467 return code+1 ;
468}
469
470
471
472////////////////////////////////////////////////////////////////////////////////
473/// Return analytical integration defined by given code
474
475double RooSimultaneous::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* /*rangeName*/) const
476{
477 // No integration scenario
478 if (code==0) {
479 return getVal(normSet) ;
480 }
481
482 // Partial integration scenarios, rangeName already encoded in 'code'
483 CacheElem* cache = static_cast<CacheElem*>(_partIntMgr.getObjByIndex(code-1)) ;
484
486 Int_t idx = _pdfProxyList.IndexOf(proxy) ;
487 return (static_cast<RooAbsReal*>(cache->_partIntList.at(idx)))->getVal(normSet) ;
488}
489
490
491
492
493
494
495////////////////////////////////////////////////////////////////////////////////
496/// Back-end for plotOn() implementation on RooSimultaneous which
497/// needs special handling because a RooSimultaneous PDF cannot
498/// project out its index category via integration. plotOn() will
499/// abort if this is requested without providing a projection dataset.
500
502{
503 // Sanity checks
504 if (plotSanityChecks(frame)) return frame ;
505
506 // Special case: if an asymmetry is requested with respect to our index
507 // category, we cannot reroute the plotting to the component pdfs. The
508 // component pdfs don't depend on the index category, so the asymmetry engine
509 // in the base class would not be able to split them by index state. Instead,
510 // we delegate directly to the base class implementation, which constructs the
511 // asymmetry from the two index-state component pdfs (see the overridden
512 // createAsymmetryComponent() and GitHub issue #14255).
513 if (auto *asymCmd = static_cast<RooCmdArg *>(cmdList.FindObject("Asymmetry"))) {
514 auto *asymCat = dynamic_cast<RooAbsCategory const *>(asymCmd->getObject(0));
515 if (asymCat && asymCat == &_indexCat.arg()) {
516
517 RooLinkedList cmdList2(cmdList);
518
519 // The base-class asymmetry-plotting engine averages the projection over
520 // the projection dataset. This is not supported for the composite data
521 // stores that back datasets with a category index, so we flatten such a
522 // projection dataset into a plain (vector-backed) copy first. Both the
523 // copy and the replacement command must outlive the plotOn() call below,
524 // because the command list only stores pointers to them.
525 std::unique_ptr<RooAbsData> flatProjData;
526 RooCmdArg newProjWData;
527 if (auto *projWData = static_cast<RooCmdArg *>(cmdList2.FindObject("ProjData"))) {
528 auto *projData = dynamic_cast<RooDataSet const *>(projWData->getObject(1));
529 if (projData && dynamic_cast<RooCompositeDataStore const *>(projData->store())) {
530 flatProjData = std::make_unique<RooDataSet>(projData->GetName(), projData->GetTitle(), *projData->get(),
531 RooFit::Import(*const_cast<RooDataSet *>(projData)));
532 const RooArgSet *projDataSet = projWData->getSet(0);
533 newProjWData = projDataSet ? RooFit::ProjWData(*projDataSet, *flatProjData)
534 : RooFit::ProjWData(*flatProjData);
535 replaceOrAdd(cmdList2, newProjWData);
536 }
537 }
538
539 return RooAbsReal::plotOn(frame, cmdList2);
540 }
541 }
542
543 // In switch mode, i.e., if the index category is neither in the projection
544 // data nor an observable of the frame, plot the selected component. This is
545 // the case for a nested RooSimultaneous when plotting a slice of the outer one.
546 if (auto *projWData = static_cast<RooCmdArg *>(cmdList.FindObject("ProjData"))) {
547 auto *projData = static_cast<RooAbsData const *>(projWData->getObject(1));
548 const RooArgSet *normVars = frame->getNormVars();
549 if (projData && !indexCatIsObservable(*projData->get()) && !(normVars && indexCatIsObservable(*normVars))) {
550 if (RooAbsPdf *pdf = getPdf(_indexCat.label())) {
551 return pdf->plotOn(frame, cmdList);
552 }
553 }
554 }
555
556 // Extract projection configuration from command list
557 RooCmdConfig pc("RooSimultaneous::plotOn(" + std::string(GetName()) + ")");
558 pc.defineString("sliceCatState","SliceCat",0,"",true) ;
559 pc.defineDouble("scaleFactor","Normalization",0,1.0) ;
560 pc.defineInt("scaleType","Normalization",0,RooAbsPdf::Relative) ;
561 pc.defineObject("sliceCatList","SliceCat",0,nullptr,true) ;
562 // This dummy is needed for plotOn to recognize the "SliceCatMany" command.
563 // It is not used directly, but the "SliceCat" commands are nested in it.
564 // Removing this dummy definition results in "ERROR: unrecognized command: SliceCatMany".
565 pc.defineObject("dummy1","SliceCatMany",0) ;
566 pc.defineSet("projSet","Project",0) ;
567 pc.defineSet("sliceSet","SliceVars",0) ;
568 pc.defineSet("projDataSet","ProjData",0) ;
569 pc.defineObject("projData","ProjData",1) ;
570 pc.defineMutex("Project","SliceVars") ;
571 pc.allowUndefined() ; // there may be commands we don't handle here
572
573 // Process and check varargs
574 pc.process(cmdList) ;
575 if (!pc.ok(true)) {
576 return frame ;
577 }
578
579 RooAbsData* projData = static_cast<RooAbsData*>(pc.getObject("projData")) ;
580 const RooArgSet* projDataSet = pc.getSet("projDataSet");
581 const RooArgSet* sliceSetTmp = pc.getSet("sliceSet") ;
582 std::unique_ptr<RooArgSet> sliceSet( sliceSetTmp ? (static_cast<RooArgSet*>(sliceSetTmp->Clone())) : nullptr );
583 const RooArgSet* projSet = pc.getSet("projSet") ;
584 double scaleFactor = pc.getDouble("scaleFactor") ;
585 ScaleType stype = (ScaleType) pc.getInt("scaleType") ;
586
587
588 // Look for category slice arguments and add them to the master slice list if found
589 const char* sliceCatState = pc.getString("sliceCatState",nullptr,true) ;
590 const RooLinkedList& sliceCatList = pc.getObjectList("sliceCatList") ;
591 if (sliceCatState) {
592
593 // Make the master slice set if it doesnt exist
594 if (!sliceSet) {
595 sliceSet = std::make_unique<RooArgSet>();
596 }
597
598 // Prepare comma separated label list for parsing
599 auto catTokens = ROOT::Split(sliceCatState, ",");
600
601 // Loop over all categories provided by (multiple) Slice() arguments
602 unsigned int tokenIndex = 0;
603 for(auto * scat : static_range_cast<RooCategory*>(sliceCatList)) {
604 const char* slabel = tokenIndex >= catTokens.size() ? nullptr : catTokens[tokenIndex++].c_str();
605
606 if (slabel) {
607 // Set the slice position to the value indicated by slabel
608 scat->setLabel(slabel) ;
609 // Add the slice category to the master slice set
610 sliceSet->add(*scat,false) ;
611 }
612 }
613 }
614
615 // Check if we have a projection dataset
616 if (!projData) {
617 coutE(InputArguments) << "RooSimultaneous::plotOn(" << GetName() << ") ERROR: must have a projection dataset for index category" << std::endl ;
618 return frame ;
619 }
620
621 // Make list of variables to be projected
622 RooArgSet projectedVars ;
623 if (sliceSet) {
624 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
625
626 // Take out the sliced variables
627 for (const auto sliceArg : *sliceSet) {
628 RooAbsArg* arg = projectedVars.find(sliceArg->GetName()) ;
629 if (arg) {
630 projectedVars.remove(*arg) ;
631 } else {
632 coutI(Plotting) << "RooAbsReal::plotOn(" << GetName() << ") slice variable "
633 << sliceArg->GetName() << " was not projected anyway" << std::endl ;
634 }
635 }
636 } else if (projSet) {
637 makeProjectionSet(frame->getPlotVar(),projSet,projectedVars,false) ;
638 } else {
639 makeProjectionSet(frame->getPlotVar(),frame->getNormVars(),projectedVars,true) ;
640 }
641
642 bool projIndex(false) ;
643
644 if (!_indexCat.arg().isDerived()) {
645 // *** Error checking for a fundamental index category ***
646 //cout << "RooSim::plotOn: index is fundamental" << std::endl ;
647
648 // Check that the provided projection dataset contains our index variable
649 if (!projData->get()->find(_indexCat.arg().GetName())) {
650 coutE(Plotting) << "RooSimultaneous::plotOn(" << GetName() << ") ERROR: Projection over index category "
651 << "requested, but projection data set doesn't contain index category" << std::endl ;
652 return frame ;
653 }
654
655 if (projectedVars.find(_indexCat.arg().GetName())) {
656 projIndex=true ;
657 }
658
659 } else {
660 // *** Error checking for a composite index category ***
661
662 // Determine if any servers of the index category are in the projectedVars
663 RooArgSet projIdxServers ;
664 bool anyServers(false) ;
665 for (const auto server : flattenedCatList()) {
666 if (projectedVars.find(server->GetName())) {
667 anyServers=true ;
668 projIdxServers.add(*server) ;
669 }
670 }
671
672 // Check that the projection dataset contains all the
673 // index category components we're projecting over
674
675 // Determine if all projected servers of the index category are in the projection dataset
676 bool allServers(true) ;
677 std::string missing;
678 for (const auto server : projIdxServers) {
679 if (!projData->get()->find(server->GetName())) {
680 allServers=false ;
681 missing = server->GetName();
682 }
683 }
684
685 if (!allServers) {
686 coutE(Plotting) << "RooSimultaneous::plotOn(" << GetName()
687 << ") ERROR: Projection dataset doesn't contain complete set of index categories to do projection."
688 << "\n\tcategory " << missing << " is missing." << std::endl ;
689 return frame ;
690 }
691
692 if (anyServers) {
693 projIndex = true ;
694 }
695 }
696
697 // Calculate relative weight fractions of components
698 std::unique_ptr<Roo1DTable> wTable( projData->table(_indexCat.arg()) );
699
700 // Clone the index category to be able to cycle through the category states for plotting without
701 // affecting the category state of our instance
702 RooArgSet idxCloneSet;
703 RooArgSet(*_indexCat).snapshot(idxCloneSet, true);
704 auto idxCatClone = static_cast<RooAbsCategoryLValue*>(idxCloneSet.find(_indexCat->GetName()) );
705 assert(idxCatClone);
706
707 // Make list of category columns to exclude from projection data
708 std::unique_ptr<RooArgSet> idxCompSliceSet( idxCatClone->getObservables(frame->getNormVars()) );
709
710 // If we don't project over the index, just do the regular plotOn
711 if (!projIndex) {
712
713 coutI(Plotting) << "RooSimultaneous::plotOn(" << GetName() << ") plot on " << frame->getPlotVar()->GetName()
714 << " represents a slice in the index category (" << _indexCat.arg().GetName() << ")" << std::endl ;
715
716 // Reduce projData: take out fitCat (component) columns and entries that don't match selected slice
717 // Construct cut string to only select projection data event that match the current slice
718
719 // Make cut string to exclude rows from projection data
720 if (sliceSet) {
721 for (auto *idxComp : static_range_cast<RooCategory *>(*idxCompSliceSet)) {
722 if (auto* slicedComponent = static_cast<const RooAbsCategory*>(sliceSet->find(*idxComp))) {
723 idxComp->setIndex(slicedComponent->getCurrentIndex(), false);
724 }
725 }
726 }
727 std::string cutString = RooFit::Detail::makeSliceCutString(*idxCompSliceSet);
728
729 // Make temporary projData without RooSim index category components
730 RooArgSet projDataVars(*projData->get()) ;
731 projDataVars.remove(*idxCompSliceSet,true,true) ;
732
733 std::unique_ptr<RooAbsData>
734 projDataTmp(projData->reduce(RooFit::SelectVars(projDataVars), RooFit::Cut(cutString.c_str())));
735
736 // Override normalization and projection dataset
737 RooCmdArg tmp1 =
738 RooFit::Normalization(scaleFactor * wTable->get(idxCatClone->getCurrentLabel()), RooAbsReal::NumEvent);
739 RooCmdArg tmp2 = RooFit::ProjWData(*projDataSet,*projDataTmp) ;
740
741 // WVE -- do not adjust normalization for asymmetry plots
742 RooLinkedList cmdList2(cmdList) ;
743 if (!cmdList.find("Asymmetry")) {
744 replaceOrAdd(cmdList2, tmp1);
745 }
746 replaceOrAdd(cmdList2, tmp2);
747
748 // Plot single component
749 RooPlot* retFrame = getPdf(idxCatClone->getCurrentLabel())->plotOn(frame,cmdList2);
750 return retFrame ;
751 }
752
753 // If we project over the index, plot using a temporary RooAddPdf
754 // using the weights from the data as coefficients
755
756 // Build the list of indexCat components that are sliced
757 idxCompSliceSet->remove(projectedVars,true,true) ;
758
759 // Make a new expression that is the weighted sum of requested components
760 RooArgList pdfCompList ;
761 RooArgList wgtCompList ;
762//RooAbsPdf* pdf ;
763 double sumWeight(0) ;
764 for(auto * proxy : static_range_cast<RooRealProxy*>(_pdfProxyList)) {
765
766 idxCatClone->setLabel(proxy->name()) ;
767
768 // Determine if this component is the current slice (if we slice)
769 bool skip(false) ;
770 for (const auto idxSliceCompArg : *idxCompSliceSet) {
771 const auto idxSliceComp = static_cast<RooAbsCategory*>(idxSliceCompArg);
772 RooAbsCategory* idxComp = static_cast<RooAbsCategory*>(idxCloneSet.find(idxSliceComp->GetName())) ;
773 if (idxComp->getCurrentIndex()!=idxSliceComp->getCurrentIndex()) {
774 skip=true ;
775 break ;
776 }
777 }
778 if (skip) continue ;
779
780 // Instantiate a RRV holding this pdfs weight
781 wgtCompList.addOwned(std::make_unique<RooRealVar>(proxy->name(),"coef",wTable->get(proxy->name())));
782 sumWeight += wTable->getFrac(proxy->name()) ;
783
784 // Add the PDF to list list
785 pdfCompList.add(proxy->arg()) ;
786 }
787
788 TString plotVarName(GetName()) ;
789 RooAddPdf plotVar{plotVarName,"weighted sum of RS components",pdfCompList,wgtCompList};
790
791 // Fix appropriate coefficient normalization in plot function
792 if (!_plotCoefNormSet.empty()) {
793 plotVar.fixAddCoefNormalization(_plotCoefNormSet) ;
794 }
795
796 std::unique_ptr<RooAbsData> projDataTmp;
797 RooArgSet projSetTmp ;
798 if (projData) {
799
800 // Construct cut string to only select projection data event that match the current slice
801 std::string cutString = RooFit::Detail::makeSliceCutString(*idxCompSliceSet);
802
803 // Make temporary projData without RooSim index category components
804 RooArgSet projDataVars(*projData->get()) ;
805 RooArgSet idxCatServers;
806 _indexCat.arg().getObservables(frame->getNormVars(), idxCatServers) ;
807
808 projDataVars.remove(idxCatServers,true,true) ;
809
810 projDataTmp = std::unique_ptr<RooAbsData>{projData->reduce(RooFit::SelectVars(projDataVars), RooFit::Cut(cutString.c_str()))};
811
812
813
814 if (projSet) {
815 projSetTmp.add(*projSet) ;
816 projSetTmp.remove(idxCatServers,true,true);
817 }
818 }
819
820
821 if (_indexCat.arg().isDerived() && !idxCompSliceSet->empty()) {
822 coutI(Plotting) << "RooSimultaneous::plotOn(" << GetName() << ") plot on " << frame->getPlotVar()->GetName()
823 << " represents a slice in index category components " << *idxCompSliceSet << std::endl ;
824
825 RooArgSet idxCompProjSet;
826 _indexCat.arg().getObservables(frame->getNormVars(), idxCompProjSet) ;
827 idxCompProjSet.remove(*idxCompSliceSet,true,true) ;
828 if (!idxCompProjSet.empty()) {
829 coutI(Plotting) << "RooSimultaneous::plotOn(" << GetName() << ") plot on " << frame->getPlotVar()->GetName()
830 << " averages with data index category components " << idxCompProjSet << std::endl ;
831 }
832 } else {
833 coutI(Plotting) << "RooSimultaneous::plotOn(" << GetName() << ") plot on " << frame->getPlotVar()->GetName()
834 << " averages with data index category (" << _indexCat.arg().GetName() << ")" << std::endl ;
835 }
836
837
838 // Override normalization and projection dataset
839 RooLinkedList cmdList2(cmdList) ;
840
841 RooCmdArg tmp1 = RooFit::Normalization(scaleFactor*sumWeight,stype) ;
842 RooCmdArg tmp2 = RooFit::ProjWData(*projDataSet,*projDataTmp) ;
843 // WVE -- do not adjust normalization for asymmetry plots
844 if (!cmdList.find("Asymmetry")) {
845 replaceOrAdd(cmdList2, tmp1);
846 }
847 replaceOrAdd(cmdList2, tmp2);
848
849 RooPlot* frame2 ;
850 if (!projSetTmp.empty()) {
851 // Plot temporary function
852 RooCmdArg tmp3 = RooFit::Project(projSetTmp) ;
853 replaceOrAdd(cmdList2, tmp3);
854 frame2 = plotVar.plotOn(frame,cmdList2) ;
855 } else {
856 // Plot temporary function
857 frame2 = plotVar.plotOn(frame,cmdList2) ;
858 }
859
860 return frame2 ;
861}
862
863
864////////////////////////////////////////////////////////////////////////////////
865/// Build the component function of an asymmetry plot (see
866/// RooAbsReal::plotAsymOn()) for a fixed state of the asymmetry category.
867///
868/// When the asymmetry is requested in our own index category, the component for
869/// a given index state is simply the corresponding pdf. We return a clone of
870/// that pdf directly instead of a RooSimultaneous with a pinned index, because
871/// a RooSimultaneous compiles its per-category observables with a category
872/// prefix. That prefix makes it incompatible with the vectorized evaluation
873/// backend that averages the asymmetry over the projection data, and would
874/// otherwise silently yield a flat (zero) asymmetry (see issue #14255). For any
875/// other asymmetry category we fall back to the generic implementation.
876
877std::unique_ptr<RooAbsReal>
879{
880 if (&asymCat == &_indexCat.arg()) {
881 const std::string &label = _indexCat.arg().lookupName(asymCatState.getCurrentIndex());
882 if (RooAbsPdf *pdf = getPdf(label)) {
883 return RooHelpers::cloneTreeWithSameParameters(static_cast<RooAbsReal const &>(*pdf));
884 }
885 }
886 return RooAbsReal::createAsymmetryComponent(asymCat, asymCatState);
887}
888
889
890////////////////////////////////////////////////////////////////////////////////
891/// Interface function used by test statistics to freeze choice of observables
892/// for interpretation of fraction coefficients. Needed here because a RooSimultaneous
893/// works like a RooAddPdf when plotted
894
895void RooSimultaneous::selectNormalization(const RooArgSet* normSet, bool /*force*/)
896{
898 if (normSet) {
899 // The index category must not be stored in the set: it is meaningless for
900 // the coefficient normalization, since it is never an observable of the
901 // component pdfs (RooAddPdf::selectNormalization() would filter it out
902 // again anyway). Worse, it is already registered as a value server via
903 // the index category proxy, and registering the same server a second
904 // time through this non-propagating set proxy corrupts the reference
905 // counts of the server's client lists when the set is cleared again.
906 RooArgSet filteredNormSet{*normSet};
907 filteredNormSet.remove(_indexCat.arg(), true, true);
908 _plotCoefNormSet.add(filteredNormSet);
909 }
910}
911
912
913////////////////////////////////////////////////////////////////////////////////
914/// Interface function used by test statistics to freeze choice of range
915/// for interpretation of fraction coefficients. Needed here because a RooSimultaneous
916/// works like a RooAddPdf when plotted
917
918void RooSimultaneous::selectNormalizationRange(const char* normRange2, bool /*force*/)
919{
920 _plotCoefNormRange = RooNameReg::ptr(normRange2) ;
921}
922
923
924
925
926////////////////////////////////////////////////////////////////////////////////
927
929 const RooArgSet* auxProto, bool verbose, bool autoBinned, const char* binnedTag) const
930{
931 const char* idxCatName = _indexCat.arg().GetName() ;
932
933 if (vars.find(idxCatName) && prototype==nullptr
934 && (auxProto==nullptr || auxProto->empty())
935 && (autoBinned || (binnedTag && strlen(binnedTag)))) {
936
937 // Return special generator config that can also do binned generation for selected states
938 return new RooSimSplitGenContext(*this,vars,verbose,autoBinned,binnedTag) ;
939
940 } else {
941
942 // Return regular generator config ;
943 return genContext(vars,prototype,auxProto,verbose) ;
944 }
945}
946
947
948
949////////////////////////////////////////////////////////////////////////////////
950/// Return specialized generator context for simultaneous p.d.f.s
951
953 const RooArgSet* auxProto, bool verbose) const
954{
955 RooArgSet allVars{vars};
956 if(prototype) allVars.add(*prototype->get());
957
958 RooArgSet catsAmongAllVars;
959 allVars.selectCommon(flattenedCatList(), catsAmongAllVars);
960
961 // Not generating index cat: we better error out because it's not clear what
962 // the user expects here. Does the user want to generate according to the
963 // currently-selected pdf? Or does the user want to generate global
964 // observable values according to the union of all category pdfs?
965 // Print an error and tell the user what to do to explicitly.
966 if(catsAmongAllVars.empty()) {
967 coutE(InputArguments) << "RooSimultaneous::generateSimGlobal(" << GetName()
968 << ") asking to generate without the index category!\n"
969 << "It's not clear what to do. you probably want to either:\n"
970 << "\n"
971 << " 1. Generate according to the currently-selected pdf.\n"
972 << " Please do this explicitly with:\n"
973 << " simpdf->getPdf(simpdf->indexCat().getCurrentLabel())->generate(vars, ...)\n"
974 << "\n"
975 << " 1. Generate global observable values according to the union of all component pdfs.\n"
976 << " For this, please use simpdf->generateSimGlobal(vars, ...)\n"
977 << std::endl;
978 return nullptr;
979 }
980
981 RooArgSet catsAmongProtoVars;
982 if(prototype) {
983 prototype->get()->selectCommon(flattenedCatList(), catsAmongProtoVars);
984
985 if(!catsAmongProtoVars.empty() && catsAmongProtoVars.size() != flattenedCatList().size()) {
986 // Abort if we have only part of the servers
987 coutE(Plotting) << "RooSimultaneous::genContext: ERROR: prototype must include either all "
988 << " components of the RooSimultaneous index category or none " << std::endl;
989 return nullptr;
990 }
991 }
992
993 return new RooSimGenContext(*this,vars,prototype,auxProto,verbose) ;
994}
995
996
997
998
999////////////////////////////////////////////////////////////////////////////////
1000
1002 const RooArgSet* nset,
1003 double scaleFactor,
1004 bool correctForBinVolume,
1005 bool showProgress) const
1006{
1007 if (RooAbsReal::fillDataHist (hist, nset, scaleFactor,
1008 correctForBinVolume, showProgress) == nullptr)
1009 return nullptr;
1010
1011 const double sum = hist->sumEntries();
1012 if (sum != 0) {
1013 for (int i=0 ; i<hist->numEntries() ; i++) {
1014 hist->set(i, hist->weight(i) / sum, 0.);
1015 }
1016 }
1017
1018 return hist;
1019}
1020
1021
1022
1023
1024////////////////////////////////////////////////////////////////////////////////
1025/// Special generator interface for generation of 'global observables' -- for RooStats tools.
1026///
1027/// \note Why one can't just use RooAbsPdf::generate()? That's because when
1028/// using the regular generate() method, a specific component pdf is selected
1029/// for each generated entry according to the index category value. However,
1030/// global observable values are independent of the current index category,
1031/// which can best be illustrated with the case where a global observable
1032/// corresponds to a nuisance parameter that is relevant for multiple channels.
1033/// So the interpretation of what is an entry in the generated dataset is very
1034/// different, hence the separate function.
1035
1037{
1038 // Generating the index category together with the global observables doesn't make any sense.
1039 RooArgSet catsAmongAllVars;
1040 whatVars.selectCommon(flattenedCatList(), catsAmongAllVars);
1041 if(!catsAmongAllVars.empty()) {
1042 coutE(InputArguments) << "RooSimultaneous::generateSimGlobal(" << GetName()
1043 << ") asking to generate global obserables at the same time as the index category!\n"
1044 << "This doesn't make any sense: global observables are generally not related to a specific channel.\n"
1045 << std::endl;
1046 return nullptr;
1047 }
1048
1049 // Make set with clone of variables (placeholder for output)
1050 RooArgSet globClone;
1051 whatVars.snapshot(globClone);
1052
1053 auto data = std::make_unique<RooDataSet>("gensimglobal","gensimglobal",whatVars);
1054
1055 for (Int_t i=0 ; i<nEvents ; i++) {
1056 for (const auto& nameIdx : indexCat()) {
1057
1058 // Get pdf associated with state from simpdf
1059 RooAbsPdf* pdftmp = getPdf(nameIdx.first);
1060
1061 RooArgSet globtmp;
1062 pdftmp->getObservables(&whatVars, globtmp) ;
1063
1064 // If there are any, generate only global variables defined by the pdf
1065 // associated with this state and transfer values to output placeholder.
1066 if (!globtmp.empty()) {
1067 globClone.assign(*std::unique_ptr<RooDataSet>{pdftmp->generate(globtmp,1)}->get(0)) ;
1068 }
1069 }
1070 data->add(globClone) ;
1071 }
1072
1073 return RooFit::makeOwningPtr(std::move(data));
1074}
1075
1076
1077/// Wraps the components of this RooSimultaneous in RooBinSamplingPdfs.
1078/// \param[in] data The dataset to be used in the eventual fit, used to figure
1079/// out the observables and whether the dataset is binned.
1080/// \param[in] precision Precision argument for all created RooBinSamplingPdfs.
1082
1083 if (precision < 0.) return;
1084
1085 RooArgSet newSamplingPdfs;
1086
1087 for (auto const &item : this->indexCat()) {
1088
1089 auto const &catName = item.first;
1090 auto &pdf = *this->getPdf(catName);
1091
1092 if (auto newSamplingPdf = RooBinSamplingPdf::create(pdf, data, precision)) {
1093 // Set the "ORIGNAME" attribute the indicate to
1094 // RooAbsArg::redirectServers() which pdf should be replaced by this
1095 // RooBinSamplingPdf in the RooSimultaneous.
1096 newSamplingPdf->setAttribute(
1097 (std::string("ORIGNAME:") + pdf.GetName()).c_str());
1098 newSamplingPdfs.addOwned(std::move(newSamplingPdf));
1099 }
1100 }
1101
1102 this->redirectServers(newSamplingPdfs, false, true);
1103 this->addOwnedComponents(std::move(newSamplingPdfs));
1104}
1105
1106
1107/// Wraps the components of this RooSimultaneous in RooBinSamplingPdfs, with a
1108/// different precision parameter for each component.
1109/// \param[in] data The dataset to be used in the eventual fit, used to figure
1110/// out the observables and whether the dataset is binned.
1111/// \param[in] precisions The map that gives the precision argument for each
1112/// component in the RooSimultaneous. The keys are the pdf names. If
1113/// there is no value for a given component, it will not use the bin
1114/// integration. Otherwise, the value has the same meaning than in
1115/// the IntegrateBins() command argument for RooAbsPdf::fitTo().
1116/// \param[in] useCategoryNames If this flag is set, the category names will be
1117/// used to look up the precision in the precisions map instead of
1118/// the pdf names.
1120 std::map<std::string, double> const& precisions,
1121 bool useCategoryNames /*=false*/) {
1122
1123 constexpr double defaultPrecision = -1.;
1124
1125 RooArgSet newSamplingPdfs;
1126
1127 for (auto const &item : this->indexCat()) {
1128
1129 auto const &catName = item.first;
1130 auto &pdf = *this->getPdf(catName);
1131 std::string pdfName = pdf.GetName();
1132
1133 auto found = precisions.find(useCategoryNames ? catName : pdfName);
1134 const double precision =
1135 found != precisions.end() ? found->second : defaultPrecision;
1136 if (precision < 0.)
1137 continue;
1138
1139 if (auto newSamplingPdf = RooBinSamplingPdf::create(pdf, data, precision)) {
1140 // Set the "ORIGNAME" attribute the indicate to
1141 // RooAbsArg::redirectServers() which pdf should be replaced by this
1142 // RooBinSamplingPdf in the RooSimultaneous.
1143 newSamplingPdf->setAttribute(
1144 (std::string("ORIGNAME:") + pdf.GetName()).c_str());
1145 newSamplingPdfs.addOwned(std::move(newSamplingPdf));
1146 }
1147 }
1148
1149 this->redirectServers(newSamplingPdfs, false, true);
1150 this->addOwnedComponents(std::move(newSamplingPdfs));
1151}
1152
1153/// Internal utility function to get a list of all category components for this
1154/// RooSimultaneous. The output contains only the index category if it is a
1155/// RooCategory, or the list of all category components if it is a
1156/// RooSuperCategory.
1158{
1159 // Note that the index category of a RooSimultaneous can only be of type
1160 // RooCategory or RooSuperCategory, because these are the only classes that
1161 // inherit from RooAbsCategoryLValue.
1162 if (auto superCat = dynamic_cast<RooSuperCategory const*>(&_indexCat.arg())) {
1163 return superCat->inputCatList();
1164 }
1165
1166 if(!_indexCatSet) {
1167 _indexCatSet = std::make_unique<RooArgSet>(_indexCat.arg());
1168 }
1169 return *_indexCatSet;
1170}
1171
1172////////////////////////////////////////////////////////////////////////////////
1173/// For the generator contexts: a component that is a RooSimultaneous whose
1174/// index category is not among the generated variables `vars` is in switch
1175/// mode, so its selected component is generated instead. Returns nullptr if
1176/// the selected state has no pdf.
1178{
1179 while (auto *sim = dynamic_cast<RooSimultaneous *>(pdf)) {
1180 if (sim->indexCatIsObservable(vars)) {
1181 break;
1182 }
1183 pdf = sim->getPdf(sim->indexCat().getCurrentLabel());
1184 }
1185 return pdf;
1186}
1187
1188////////////////////////////////////////////////////////////////////////////////
1189/// Check if the index category is in the normalization set `nset`, matching by
1190/// name. For a RooSuperCategory index, it's enough if all its input categories
1191/// are in `nset`, which is the case for the variables of a dataset.
1193{
1194 if (!nset) {
1195 return false;
1196 }
1197 if (nset->contains(_indexCat.arg())) {
1198 return true;
1199 }
1200 RooArgSet const &cats = flattenedCatList();
1201 return std::all_of(cats.begin(), cats.end(), [&](RooAbsArg const *cat) { return nset->contains(*cat); });
1202}
1203
1204////////////////////////////////////////////////////////////////////////////////
1205/// Check if the index category is among the variables `vars`, matching by
1206/// name. For a RooSuperCategory index, any of its input categories counts.
1207///
1208/// If the index category is not among the observables of a fit, this
1209/// RooSimultaneous does not split the data into channels: it acts as a
1210/// "switch" that evaluates to the component selected by the current index
1211/// state, analogous to RooMultiPdf. Fitting infrastructure uses this check to
1212/// decide between the two modes.
1214{
1215 RooArgSet catsAmongVars;
1216 vars.selectCommon(flattenedCatList(), catsAmongVars);
1217 return !catsAmongVars.empty();
1218}
1219
1220namespace {
1221
1222void markObs(RooAbsArg *arg, std::string const &prefix, RooArgSet const &normSet)
1223{
1224 for (RooAbsArg *server : arg->servers()) {
1225 if (server->isFundamental() && normSet.find(*server)) {
1226 markObs(server, prefix, normSet);
1227 server->setAttribute("__obs__");
1228 } else if (!server->isFundamental()) {
1229 markObs(server, prefix, normSet);
1230 }
1231 }
1232}
1233
1234void prefixArgs(RooAbsArg *arg, std::string const &prefix, RooArgSet const &normSet)
1235{
1236 if (!arg->getStringAttribute("__prefix__")) {
1237 arg->SetName((prefix + arg->GetName()).c_str());
1238 arg->setStringAttribute("__prefix__", prefix.c_str());
1239 }
1240 for (RooAbsArg *server : arg->servers()) {
1241 if (server->isFundamental() && normSet.find(*server)) {
1242 prefixArgs(server, prefix, normSet);
1243 } else if (!server->isFundamental()) {
1244 prefixArgs(server, prefix, normSet);
1245 }
1246 }
1247}
1248
1249} // namespace
1250
1251std::unique_ptr<RooAbsArg>
1253{
1254 if (!indexCatIsObservable(normSet)) {
1255 // The index category is not an observable here: the RooSimultaneous
1256 // acts as a plain "switch" that evaluates to the component selected by
1257 // the current index state, analogous to RooMultiPdf. The channel
1258 // observables are then the same as the ones of this pdf, so the
1259 // channel-splitting compilation below (which renames the per-channel
1260 // observables so they can be filled from split datasets) must not be
1261 // used. Compile like an ordinary self-normalized pdf instead.
1262 return RooAbsPdf::compileForNormSet(normSet, ctx);
1263 }
1264
1265 std::unique_ptr<RooSimultaneous> newSimPdf{static_cast<RooSimultaneous *>(this->Clone())};
1266
1267 const char *rangeName = this->getStringAttribute("RangeName");
1268 bool splitRange = this->getAttribute("SplitRange");
1269
1270 RooArgSet newPdfs;
1271 std::vector<std::string> catNames;
1272
1273 for (auto *proxy : static_range_cast<RooRealProxy *>(newSimPdf->_pdfProxyList)) {
1274 catNames.emplace_back(proxy->GetName());
1275 std::string const &catName = catNames.back();
1276 const std::string prefix = "_" + catName + "_";
1277
1278 const std::string origname = proxy->arg().GetName();
1279
1280 auto pdfClone = RooHelpers::cloneTreeWithSameParameters(static_cast<RooAbsPdf const &>(proxy->arg()), &normSet);
1281
1282 markObs(pdfClone.get(), prefix, normSet);
1283
1284 std::unique_ptr<RooArgSet> pdfNormSet{
1285 std::unique_ptr<RooArgSet>(pdfClone->getVariables())->selectByAttrib("__obs__", true)};
1286 std::unique_ptr<RooArgSet> condVarSet{
1287 std::unique_ptr<RooArgSet>(pdfClone->getVariables())->selectByAttrib("__conditional__", true)};
1288
1289 pdfNormSet->remove(*condVarSet, true, true);
1290
1291 if (rangeName) {
1292 pdfClone->setNormRange(RooHelpers::getRangeNameForSimComponent(rangeName, splitRange, catName).c_str());
1293 }
1294
1295 RooFit::Detail::CompileContext pdfContext{*pdfNormSet};
1296 pdfContext.setLikelihoodMode(ctx.likelihoodMode());
1297 auto *pdfFinal = pdfContext.compile(*pdfClone, *newSimPdf, *pdfNormSet);
1298
1299 // The observables of the compiled computation graph are clones of the
1300 // ones in `pdfNormSet`, and only the clones are prefixed below. That's
1301 // why we have to translate the normalization set to the compiled args
1302 // before using it any further.
1303 RooArgSet compiledPdfNormSet = pdfContext.mapToCompiled(*pdfNormSet);
1304
1305 // We can only prefix the observables after everything related the
1306 // compiling of the compute graph for the normalization set is done. This
1307 // is because of a subtlety in conditional RooProdPdfs, which stores the
1308 // normalization sets for the individual pdfs in RooArgSets that are
1309 // disconnected from the computation graph, so we have no control over
1310 // them. An alternative would be to use recursive server re-direction,
1311 // but this has more performance overhead.
1312 prefixArgs(pdfFinal, prefix, normSet);
1313
1314 // This has to be done with the observables of the compiled graph, which
1315 // are prefixed at this point. Otherwise, the RooAddPdf components don't
1316 // find their coefficient normalization observables and silently keep the
1317 // default "automatic" interpretation of the coefficients, which results
1318 // in wrong yields in ranged fits (GitHub issue #23444).
1319 pdfFinal->fixAddCoefNormalization(compiledPdfNormSet, false);
1320
1321 pdfClone->SetName((std::string("_") + pdfClone->GetName()).c_str());
1322 pdfFinal->addOwnedComponents(std::move(pdfClone));
1323
1324 pdfFinal->setAttribute(("ORIGNAME:" + origname).c_str());
1325 newPdfs.add(*pdfFinal);
1326
1327 // We will remove the old pdf server because we will fill the new ones by
1328 // hand via the creation of new proxies.
1329 newSimPdf->removeServer(const_cast<RooAbsReal &>(proxy->arg()), true);
1330 }
1331
1332 // Replace pdfs with compiled pdfs. Don't use RooAbsArg::redirectServers()
1333 // here, because it doesn't support replacing two servers with the same name
1334 // (it can happen in a RooSimultaneous that two pdfs have the same name).
1335
1336 // First delete old proxies (we have already removed the servers before).
1337 newSimPdf->_pdfProxyList.Delete();
1338
1339 // Recreate the _pdfProxyList with the compiled pdfs
1340 for (std::size_t i = 0; i < newPdfs.size(); ++i) {
1341 const char *label = catNames[i].c_str();
1342 newSimPdf->_pdfProxyList.Add(
1343 new RooRealProxy(label, label, newSimPdf.get(), *static_cast<RooAbsReal *>(newPdfs[i])));
1344 }
1345
1346 ctx.compileServers(*newSimPdf, normSet); // to trigger compiling also the index category
1347
1348 return newSimPdf;
1349}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define coutE(a)
RooTemplateProxy< RooAbsReal > RooRealProxy
Compatibility typedef replacing the old RooRealProxy class.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition TGX11.cxx:142
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
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.
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
void SetName(const char *name) override
Set the name of the TNamed.
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
bool redirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool isRecursionStep=false)
Replace all direct servers of this object with the new servers in newServerList.
virtual bool isDerived() const
Does value or shape of this arg depend on any other arg?
Definition RooAbsArg.h:98
bool getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Definition RooAbsArg.h:89
Abstract base class for objects that represent a discrete value that can be set from the outside,...
A space to attach TBranches.
virtual value_type getCurrentIndex() const
Return index number of current state.
const std::string & lookupName(value_type index) const
Get the name corresponding to the given index.
std::map< std::string, value_type >::const_iterator begin() const
Iterator for category state names. Points to pairs of index and name.
std::size_t size() const
Number of states defined.
bool contains(const char *name) const
Check if collection contains an argument with a specific name.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
const_iterator end() const
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...
Storage_t::size_type size() const
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
const_iterator begin() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:55
virtual const RooArgSet * get() const
Definition RooAbsData.h:99
virtual Roo1DTable * table(const RooArgSet &catSet, const char *cuts="", const char *opts="") const
Construct table for product of categories in catSet.
RooFit::OwningPtr< RooAbsData > reduce(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a reduced copy of this dataset.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract base class for generator contexts of RooAbsPdf objects.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:32
std::unique_ptr< RooAbsArg > compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const override
RooArgSet const * _normSet
! Normalization set with for above integral
Definition RooAbsPdf.h:314
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:116
bool canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:214
@ CanBeExtended
Definition RooAbsPdf.h:208
@ MustBeExtended
Definition RooAbsPdf.h:208
@ CanNotBeExtended
Definition RooAbsPdf.h:208
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:49
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:63
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.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:107
bool plotSanityChecks(RooPlot *frame) const
Utility function for plotOn(), perform general sanity check on frame to ensure safe plotting operatio...
virtual std::unique_ptr< RooAbsReal > createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const
Build the component function of an asymmetry plot (see plotAsymOn()) that corresponds to a fixed stat...
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const
Plot (project) PDF on specified frame.
void makeProjectionSet(const RooAbsArg *plotVar, const RooArgSet *allVars, RooArgSet &projectedVars, bool silent) const
Utility function for plotOn() that constructs the set of observables to project when plotting ourselv...
Efficient implementation of a sum of PDFs of the form.
Definition RooAddPdf.h:32
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:24
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:159
RooArgSet * selectCommon(const RooAbsCollection &refColl) const
Use RooAbsCollection::selecCommon(), but return as RooArgSet.
Definition RooArgSet.h:154
static std::unique_ptr< RooAbsPdf > create(RooAbsPdf &pdf, RooAbsData const &data, double precision)
Creates a wrapping RooBinSamplingPdf if appropriate.
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
T * getObjByIndex(Int_t index) const
Retrieve payload object by slot index.
Int_t lastIndex() const
Return index of slot used in last get or set operation.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter function without integration set.
Named container for two doubles, two integers two object points and three string pointers that can be...
Definition RooCmdArg.h:26
Configurable parser for RooCmdArg named arguments.
void defineMutex(const char *head, Args_t &&... tail)
Define arguments where any pair is mutually exclusive.
bool process(const RooCmdArg &arg)
Process given RooCmdArg.
double getDouble(const char *name, double defaultValue=0.0) const
Return double property registered with name 'name'.
bool defineDouble(const char *name, const char *argName, int doubleNum, double defValue=0.0)
Define double property name 'name' mapped to double in slot 'doubleNum' in RooCmdArg with name argNam...
RooArgSet * getSet(const char *name, RooArgSet *set=nullptr) const
Return RooArgSet property registered with name 'name'.
bool defineSet(const char *name, const char *argName, int setNum, const RooArgSet *set=nullptr)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
bool ok(bool verbose) const
Return true of parsing was successful.
bool defineObject(const char *name, const char *argName, int setNum, const TObject *obj=nullptr, bool isArray=false)
Define TObject property name 'name' mapped to object in slot 'setNum' in RooCmdArg with name argName ...
const char * getString(const char *name, const char *defaultValue="", bool convEmptyToNull=false) const
Return string property registered with name 'name'.
bool defineString(const char *name, const char *argName, int stringNum, const char *defValue="", bool appendMode=false)
Define double property name 'name' mapped to double in slot 'stringNum' in RooCmdArg with name argNam...
const RooLinkedList & getObjectList(const char *name) const
Return list of objects registered with name 'name'.
bool defineInt(const char *name, const char *argName, int intNum, int 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.
int getInt(const char *name, int defaultValue=0) const
Return integer property registered with name 'name'.
TObject * getObject(const char *name, TObject *obj=nullptr) const
Return TObject property registered with name 'name'.
void removeAll() override
Remove all argument inset using remove(const RooAbsArg&).
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
Combines several disjunct datasets into one.
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
double weight(std::size_t i) const
Return weight of i-th bin.
void set(std::size_t binNumber, double weight, double wgtErr)
Set bin content of bin that was last loaded with get(std::size_t).
double sumEntries() const override
Sum the weights of all bins.
Container class to hold unbinned data.
Definition RooDataSet.h:32
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
void compileServers(RooAbsArg &arg, RooArgSet const &normSet)
Collection class for internal use, storing a collection of RooAbsArg pointers in a doubly linked list...
bool Replace(const TObject *oldArg, const TObject *newArg)
Replace object 'oldArg' in collection with new object 'newArg'.
TObject * find(const char *name) const
Return pointer to object with given name in collection.
virtual void Add(TObject *arg)
TObject * FindObject(const char *name) const override
Return pointer to object with given name.
static const TNamed * ptr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
const RooArgSet * getNormVars() const
Definition RooPlot.h:146
RooAbsRealLValue * getPlotVar() const
Definition RooPlot.h:137
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
void selectNormalization(const RooArgSet *depSet=nullptr, bool force=false) override
Interface function used by test statistics to freeze choice of observables for interpretation of frac...
double evaluate() const override
Return the current value: the value of the PDF associated with the current index category state.
Int_t _numPdf
Number of registered PDFs.
void selectNormalizationRange(const char *rangeName=nullptr, bool force=false) override
Interface function used by test statistics to freeze choice of range for interpretation of fraction c...
TList _pdfProxyList
List of PDF proxies (named after applicable category state)
std::unique_ptr< RooAbsReal > createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const override
Build the component function of an asymmetry plot (see RooAbsReal::plotAsymOn()) for a fixed state of...
bool indexCatIsInNormSet(RooArgSet const *nset) const
Check if the index category is in the normalization set nset, matching by name.
std::unique_ptr< RooAbsReal > createExpectedEventsFunc(const RooArgSet *nset) const override
Expected events function for the new evaluation backends: the sum over all components,...
RooObjCacheManager _partIntMgr
! Component normalization manager
~RooSimultaneous() override
Destructor.
RooFit::OwningPtr< RooDataSet > generateSimGlobal(const RooArgSet &whatVars, Int_t nEvents) override
Special generator interface for generation of 'global observables' – for RooStats tools.
RooArgSet const & flattenedCatList() const
Internal utility function to get a list of all category components for this RooSimultaneous.
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
ExtendMode extendMode() const override
Examine the pdf components and check if one of them can be extended or must be extended.
RooCategoryProxy _indexCat
Index category.
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Return analytical integration defined by given code.
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Forward determination of analytical integration capabilities to component p.d.f.s A unique code is as...
double expectedEvents(const RooArgSet *nset) const override
Return the number of expected events: If the index is in nset, then return the sum of the expected ev...
friend class RooSimGenContext
RooAbsGenContext * autoGenContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false, bool autoBinned=true, const char *binnedTag="") const override
RooAbsPdf * getPdf(RooStringView catName) const
Return the p.d.f associated with the given index category name.
std::unique_ptr< RooAbsArg > compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext &ctx) const override
void wrapPdfsInBinSamplingPdfs(RooAbsData const &data, double precision)
Wraps the components of this RooSimultaneous in RooBinSamplingPdfs.
static RooAbsPdf * resolveSwitchMode(RooAbsPdf *pdf, RooArgSet const &vars)
For the generator contexts: a component that is a RooSimultaneous whose index category is not among t...
bool indexCatIsObservable(RooArgSet const &vars) const
Check if the index category is among the variables vars, matching by name.
const TNamed * _plotCoefNormRange
RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const override
Return specialized generator context for simultaneous p.d.f.s.
std::unique_ptr< RooArgSet > _indexCatSet
! Index category wrapped in a RooArgSet if needed internally
bool addPdf(const RooAbsPdf &pdf, const char *catLabel)
Associate given PDF with index category state label 'catLabel'.
RooSetProxy _plotCoefNormSet
const RooAbsCategoryLValue & indexCat() const
friend class RooSimSplitGenContext
virtual RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, double scaleFactor, bool correctForBinVolume=false, bool showProgress=false) const
The RooStringView is a wrapper around a C-style string that can also be constructed from a std::strin...
const char * c_str() const
Joins several RooAbsCategoryLValue objects into a single category.
const char * label() const
Get the label of the current category state. This function only makes sense for category proxies.
const T & arg() const
Return reference to object held in proxy.
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:708
void Add(TObject *obj) override
Definition TList.h:81
void Delete(Option_t *option="") override
Remove all objects from the list AND delete all heap based objects.
Definition TList.cxx:600
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
Mother of all ROOT objects.
Definition TObject.h:42
virtual const char * GetName() const
Returns name of object.
Definition TObject.cxx:461
virtual Int_t IndexOf(const TObject *obj) const
Return index of object in collection.
Basic string class.
Definition TString.h:137
RooCmdArg SelectVars(const RooArgSet &vars)
RooCmdArg Import(const char *state, TH1 &histo)
RooCmdArg ProjWData(const RooAbsData &projData, bool binData=false)
RooCmdArg Project(const RooArgSet &projSet)
RooCmdArg Normalization(double scaleFactor)
RooCmdArg Cut(const char *cutSpec)
std::vector< std::string > Split(std::string_view str, std::string_view delims, bool skipEmpty=false)
Splits a string at each character in delims.
std::string makeSliceCutString(RooArgSet const &sliceDataSet)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:73
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
std::unique_ptr< T > cloneTreeWithSameParameters(T const &arg, RooArgSet const *observables=nullptr)
Clone RooAbsArg object and reattach to original parameters.
std::string getRangeNameForSimComponent(std::string const &rangeName, bool splitRange, std::string const &catName)
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2335