diff --git a/roofit/roofitcore/inc/RooAbsReal.h b/roofit/roofitcore/inc/RooAbsReal.h index 18ca28e276d5d..80d1053ed1bf5 100644 --- a/roofit/roofitcore/inc/RooAbsReal.h +++ b/roofit/roofitcore/inc/RooAbsReal.h @@ -496,6 +496,9 @@ class RooAbsReal : public RooAbsArg { virtual RooPlot *plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asymCat, PlotOpt o) const; + virtual std::unique_ptr + createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const; + bool matchArgsByName(const RooArgSet &allArgs, RooArgSet &matchedArgs, const TList &nameList) const; bool redirectServersHook(const RooAbsCollection & newServerList, bool mustReplaceAll, diff --git a/roofit/roofitcore/inc/RooSimultaneous.h b/roofit/roofitcore/inc/RooSimultaneous.h index f5f48383412a6..e9ce1e295800a 100644 --- a/roofit/roofitcore/inc/RooSimultaneous.h +++ b/roofit/roofitcore/inc/RooSimultaneous.h @@ -113,6 +113,9 @@ class RooSimultaneous : public RooAbsPdf { void selectNormalization(const RooArgSet* depSet=nullptr, bool force=false) override ; void selectNormalizationRange(const char* rangeName=nullptr, bool force=false) override ; + std::unique_ptr + createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const override; + RooArgSet const& flattenedCatList() const; mutable RooSetProxy _plotCoefNormSet ; diff --git a/roofit/roofitcore/src/RooAbsOptTestStatistic.cxx b/roofit/roofitcore/src/RooAbsOptTestStatistic.cxx index b9e32cff6337f..595a439638624 100644 --- a/roofit/roofitcore/src/RooAbsOptTestStatistic.cxx +++ b/roofit/roofitcore/src/RooAbsOptTestStatistic.cxx @@ -226,7 +226,6 @@ void RooAbsOptTestStatistic::initSlave(RooAbsReal& real, RooAbsData& indata, con // Copy data and strip entries lost by adjusted fit range, _dataClone ranges will be copied from realDepSet ranges if (rangeName && strlen(rangeName)) { _dataClone = std::unique_ptr{indata.reduce(RooFit::SelectVars(*_funcObsSet),RooFit::CutRange(rangeName))}.release(); - // std::cout << "RooAbsOptTestStatistic: reducing dataset to fit in range named " << rangeName << " resulting dataset has " << _dataClone->sumEntries() << " events" << std::endl ; } else { _dataClone = static_cast(indata.Clone()) ; } @@ -432,9 +431,6 @@ void RooAbsOptTestStatistic::constOptimizeTestStatistic(ConstOpCode opcode, bool hasWarned = true; } - // std::cout << "ROATS::constOpt(" << GetName() << ") funcClone structure dump BEFORE const-opt" << std::endl ; - // _funcClone->Print("t") ; - RooAbsTestStatistic::constOptimizeTestStatistic(opcode,doAlsoTrackingOpt); if (operMode()!=Slave) return ; @@ -489,8 +485,6 @@ void RooAbsOptTestStatistic::constOptimizeTestStatistic(ConstOpCode opcode, bool break ; } -// std::cout << "ROATS::constOpt(" << GetName() << ") funcClone structure dump AFTER const-opt" << std::endl ; -// _funcClone->Print("t") ; } @@ -506,8 +500,6 @@ void RooAbsOptTestStatistic::constOptimizeTestStatistic(ConstOpCode opcode, bool void RooAbsOptTestStatistic::optimizeCaching() { -// std::cout << "RooAbsOptTestStatistic::optimizeCaching(" << GetName() << "," << this << ")" << std::endl ; - // Trigger create of all object caches now in nodes that have deferred object creation // so that cache contents can be processed immediately _funcClone->getVal(_normSet) ; @@ -657,13 +649,9 @@ bool RooAbsOptTestStatistic::setDataSlave(RooAbsData& indata, bool cloneData, bo { if (operMode()==SimMaster) { - //cout << "ROATS::setDataSlave() ERROR this is SimMaster _funcClone = " << _funcClone << std::endl ; return false ; } - //cout << "ROATS::setDataSlave() new dataset size = " << indata.numEntries() << std::endl ; - //indata.Print("v") ; - // If the current dataset is owned, transfer the ownership to unique pointer // that will get out of scope at the end of this function. We can't delete it diff --git a/roofit/roofitcore/src/RooAbsReal.cxx b/roofit/roofitcore/src/RooAbsReal.cxx index 72e354575bf50..979a079c9c715 100644 --- a/roofit/roofitcore/src/RooAbsReal.cxx +++ b/roofit/roofitcore/src/RooAbsReal.cxx @@ -2203,6 +2203,12 @@ RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asym } } + // The asymmetry category itself defines the two sides of the asymmetry, so it + // must not be treated as a variable to be averaged over the projection data. + if (RooAbsArg *asymCatInProjData = projDataVars.find(asymCat.GetName())) { + projDataVars.remove(*asymCatInProjData) ; + } + // Must depend on asymCat if (!dependsOn(asymCat)) { coutE(Plotting) << "RooAbsReal::plotAsymOn(" << GetName() @@ -2273,19 +2279,16 @@ RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asym } - // Customize two copies of projection with fixed negative and positive asymmetry + // Build two copies of the function with the asymmetry category fixed to its + // negative and positive state. By default these are copies with the category + // pinned via a RooCustomizer, but subclasses (RooSimultaneous) can provide a + // more suitable construction. std::unique_ptr asymPos{static_cast(asymCat.Clone("asym_pos"))}; std::unique_ptr asymNeg{static_cast(asymCat.Clone("asym_neg"))}; asymPos->setIndex(1) ; asymNeg->setIndex(-1) ; - RooCustomizer custPos{*this,"pos"}; - RooCustomizer custNeg{*this,"neg"}; - //custPos->setOwning(true) ; - //custNeg->setOwning(true) ; - custPos.replaceArg(asymCat,*asymPos) ; - custNeg.replaceArg(asymCat,*asymNeg) ; - std::unique_ptr funcPos{static_cast(custPos.build())}; - std::unique_ptr funcNeg{static_cast(custNeg.build())}; + std::unique_ptr funcPos = createAsymmetryComponent(asymCat, *asymPos); + std::unique_ptr funcNeg = createAsymmetryComponent(asymCat, *asymNeg); // Create projection integral RooArgSet *posProjCompList; @@ -2293,8 +2296,12 @@ RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asym // Add projDataVars to normalized dependents of projection // This is needed only for asymmetries (why?) - RooArgSet depPos(*plotVar,*asymPos) ; - RooArgSet depNeg(*plotVar,*asymNeg) ; + RooArgSet depPos(*plotVar) ; + RooArgSet depNeg(*plotVar) ; + // Keep the fixed asymmetry category in the normalization set only if the + // component function actually depends on it (i.e. it was pinned in place). + if (funcPos->dependsOn(*asymPos)) depPos.add(*asymPos) ; + if (funcNeg->dependsOn(*asymNeg)) depNeg.add(*asymNeg) ; depPos.add(projDataVars) ; depNeg.add(projDataVars) ; @@ -2413,6 +2420,21 @@ RooPlot* RooAbsReal::plotAsymOn(RooPlot *frame, const RooAbsCategoryLValue& asym } +//////////////////////////////////////////////////////////////////////////////// +/// Build the component function of an asymmetry plot (see plotAsymOn()) that +/// corresponds to a fixed state of the asymmetry category. The default +/// implementation returns a copy of this function with the asymmetry category +/// pinned to the requested state via a RooCustomizer. + +std::unique_ptr +RooAbsReal::createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const +{ + RooCustomizer cust{*this, asymCatState.GetName()}; + cust.replaceArg(asymCat, asymCatState); + return std::unique_ptr{static_cast(cust.build())}; +} + + //////////////////////////////////////////////////////////////////////////////// /// \brief Propagates parameter uncertainties to an uncertainty estimate for this RooAbsReal. diff --git a/roofit/roofitcore/src/RooChangeTracker.cxx b/roofit/roofitcore/src/RooChangeTracker.cxx index 386d7e9e04638..9c3b418f2af23 100644 --- a/roofit/roofitcore/src/RooChangeTracker.cxx +++ b/roofit/roofitcore/src/RooChangeTracker.cxx @@ -113,12 +113,9 @@ bool RooChangeTracker::hasChanged(bool clearState) if (clearState) { // Clear dirty flag by calling getVal() - //cout << "RooChangeTracker(" << GetName() << ") clearing isValueDirty" << std::endl ; clearValueDirty() ; } - //cout << "RooChangeTracker(" << GetName() << ") isValueDirty = true, returning true" << std::endl ; - return true ; } @@ -131,7 +128,6 @@ bool RooChangeTracker::hasChanged(bool clearState) for (unsigned int i=0; i < _realSet.size(); ++i) { auto real = static_cast(_realSet.at(i)); if (real->getVal() != _realRef[i]) { - // std::cout << "RooChangeTracker(" << this << "," << GetName() << ") value of " << real->GetName() << " has changed from " << _realRef[i] << " to " << real->getVal() << " clearState = " << (clearState?"T":"F") << std::endl ; valuesChanged = true ; _realRef[i] = real->getVal() ; } @@ -140,7 +136,6 @@ bool RooChangeTracker::hasChanged(bool clearState) for (unsigned int i=0; i < _catSet.size(); ++i) { auto cat = static_cast(_catSet.at(i)); if (cat->getCurrentIndex() != _catRef[i]) { - // std::cout << "RooChangeTracker(" << this << "," << GetName() << ") value of " << cat->GetName() << " has changed from " << _catRef[i-1] << " to " << cat->getIndex() << std::endl ; valuesChanged = true ; _catRef[i] = cat->getCurrentIndex() ; } @@ -154,8 +149,6 @@ bool RooChangeTracker::hasChanged(bool clearState) _init = true ; } - // std::cout << "RooChangeTracker(" << GetName() << ") returning " << (valuesChanged?"T":"F") << std::endl ; - return valuesChanged ; } else { diff --git a/roofit/roofitcore/src/RooConvIntegrandBinding.cxx b/roofit/roofitcore/src/RooConvIntegrandBinding.cxx index f1c3f66da6ba7..44d5e9b7180fc 100644 --- a/roofit/roofitcore/src/RooConvIntegrandBinding.cxx +++ b/roofit/roofitcore/src/RooConvIntegrandBinding.cxx @@ -113,7 +113,6 @@ void RooConvIntegrandBinding::loadValues(const double xvector[], bool clipInvali if (clipInvalid && !_vars[index]->isValidReal(xvector[index])) { _xvecValid = false ; } else { - //cout << "RooConvBasBinding::loadValues[" << index << "] loading value " << xvector[index] << std::endl ; _vars[index]->setVal(xvector[index]); } } @@ -131,7 +130,6 @@ double RooConvIntegrandBinding::operator()(const double xvector[]) const // First evaluate function at x' loadValues(xvector); if (!_xvecValid) return 0 ; - //cout << "RooConvIntegrandBinding::operator(): evaluating f(x') at x' = " << xvector[0] << std::endl ; double f_xp = _func->getVal(_nset) ; // Next evaluate model at x-x' @@ -140,12 +138,6 @@ double RooConvIntegrandBinding::operator()(const double xvector[]) const if (!_xvecValid) return 0 ; double g_xmxp = _model->getVal(_nset) ; - //cout << "RooConvIntegrandBinding::operator(): evaluating g(x-x') at x-x' = " << _vars[0]->getVal() << " = " << g_xmxp << std::endl ; - //cout << "RooConvIntegrandBinding::operator(): return value = " << f_xp << " * " << g_xmxp << " = " << f_xp*g_xmxp << std::endl ; - - //cout << "_vars[0] = " << _vars[0]->getVal() << " _vars[1] = " << _vars[1]->getVal() << std::endl ; - //cout << "_xvec[0] = " << xvector[0] << " _xvec[1] = " << xvector[1] << std::endl ; - return f_xp*g_xmxp ; } diff --git a/roofit/roofitcore/src/RooHistFunc.cxx b/roofit/roofitcore/src/RooHistFunc.cxx index 043df4cd492e4..1ffac2b904e9d 100644 --- a/roofit/roofitcore/src/RooHistFunc.cxx +++ b/roofit/roofitcore/src/RooHistFunc.cxx @@ -348,9 +348,6 @@ std::list* RooHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo, } } - // std::cout << "RooHistFunc::bb(" << GetName() << ") histObs = " << _histObsList << std::endl ; - // std::cout << "RooHistFunc::bb(" << GetName() << ") pdfObs = " << _depList << std::endl ; - RooAbsRealLValue* transform = nullptr; if (!hobs) { @@ -377,9 +374,6 @@ std::list* RooHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo, } - // std::cout << "hobs = " << hobs->GetName() << std::endl ; - // std::cout << "transform = " << (transform?transform->GetName():"") << std::endl ; - // Check that observable is in dataset, if not no hint is generated RooAbsArg* xtmp = _dataHist->get()->find(hobs->GetName()) ; if (!xtmp) { @@ -409,7 +403,6 @@ std::list* RooHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo, double boundary = boundaries[i] ; if (transform) { transform->setVal(boundary) ; - //cout << "transform bound " << boundary << " using " << transform->GetName() << " result " << obs.getVal() << std::endl ; hint->push_back(obs.getVal()) ; } else { hint->push_back(boundary) ; diff --git a/roofit/roofitcore/src/RooLinkedList.cxx b/roofit/roofitcore/src/RooLinkedList.cxx index a16a26044c9d2..9765489a5a0f9 100644 --- a/roofit/roofitcore/src/RooLinkedList.cxx +++ b/roofit/roofitcore/src/RooLinkedList.cxx @@ -55,7 +55,6 @@ namespace RooLinkedListImplDetails { _sz(sz), _free(capacity()), _chunk(new RooLinkedListElem[_free]), _freelist(_chunk) { - //cout << "RLLID::Chunk ctor(" << this << ") of size " << _free << " list elements" << std::endl ; // initialise free list for (Int_t i = 0; i < _free; ++i) _chunk[i]._next = (i + 1 < _free) ? &_chunk[i + 1] : nullptr; @@ -285,7 +284,6 @@ RooLinkedList::RooLinkedList(const RooLinkedList& other) : } //////////////////////////////////////////////////////////////////////////////// -/// std::cout << "RooLinkedList::createElem(" << this << ") obj = " << obj << " elem = " << elem << std::endl ; RooLinkedListElem* RooLinkedList::createElement(TObject* obj, RooLinkedListElem* elem) { @@ -426,7 +424,6 @@ void RooLinkedList::Add(TObject* arg, Int_t refCount) } if (_htableName){ - //cout << "storing link " << _last << " with hash arg " << arg << std::endl ; _htableName->insert({arg->GetName(), arg}); _htableLink->insert({arg, reinterpret_cast(_last)}); } @@ -609,7 +606,6 @@ TObject* RooLinkedList::find(const char* name) const if (_useNptr) { // See if it might have been renamed const TNamed* nptr= RooNameReg::known(name); - //cout << "RooLinkedList::find: possibly renamed '" << name << "', kRenamedArg=" << (nptr&&nptr->TestBit(RooNameReg::kRenamedArg)) << std::endl; if (nptr && nptr->TestBit(RooNameReg::kRenamedArg)) { RooLinkedListElem* ptr = _first ; while(ptr) { @@ -622,7 +618,6 @@ TObject* RooLinkedList::find(const char* name) const } return nullptr ; } - //cout << "RooLinkedList::find: possibly renamed '" << name << "'" << std::endl; } RooLinkedListElem* ptr = _first ; @@ -661,7 +656,6 @@ RooAbsArg* RooLinkedList::findArg(const RooAbsArg* arg) const if (_htableName) { RooAbsArg* a = const_cast(static_cast((*_htableName)[arg->GetName()])); if (a) return a; - //cout << "RooLinkedList::findArg: possibly renamed '" << arg->GetName() << "', kRenamedArg=" << arg->namePtr()->TestBit(RooNameReg::kRenamedArg) << std::endl; // See if it might have been renamed if (!arg->namePtr()->TestBit(RooNameReg::kRenamedArg)) return nullptr; } diff --git a/roofit/roofitcore/src/RooParamBinning.cxx b/roofit/roofitcore/src/RooParamBinning.cxx index 709d41dd317dd..482c14904e732 100644 --- a/roofit/roofitcore/src/RooParamBinning.cxx +++ b/roofit/roofitcore/src/RooParamBinning.cxx @@ -73,20 +73,16 @@ RooParamBinning::~RooParamBinning() //////////////////////////////////////////////////////////////////////////////// /// Copy constructor -/// std::cout << "RooParamBinning::cctor(" << this << ") orig = " << &other << std::endl ; RooParamBinning::RooParamBinning(const RooParamBinning &other, const char *name) : RooAbsBinning(name) { if (other._lp) { -// std::cout << "RooParamBinning::cctor(this = " << this << ") taking addresses from orig ListProxy" << std::endl ; _xlo = static_cast(other._lp->at(0)) ; _xhi = static_cast(other._lp->at(1)) ; } else { -// std::cout << "RooParamBinning::cctor(this = " << this << ") taking addresses from orig pointers " << other._xlo << " " << other._xhi << std::endl ; - _xlo = other._xlo ; _xhi = other._xhi ; } @@ -94,7 +90,6 @@ RooParamBinning::RooParamBinning(const RooParamBinning &other, const char *name) _nbins = other._nbins ; _lp = nullptr ; - //cout << "RooParamBinning::cctor(this = " << this << " xlo = " << &_xlo << " xhi = " << &_xhi << " _lp = " << _lp << " owner = " << _owner << ")" << std::endl ; } @@ -110,14 +105,11 @@ void RooParamBinning::insertHook(RooAbsRealLValue& owner) const _owner = &owner ; // If list proxy already exists update pointers from proxy -// std::cout << "RooParamBinning::insertHook(" << this << "," << GetName() << ") _lp at beginning = " << _lp << std::endl ; if (_lp) { -// std::cout << "updating raw pointers from list proxy contents" << std::endl ; _xlo = xlo() ; _xhi = xhi() ; delete _lp ; } -// std::cout << "_xlo = " << _xlo << " _xhi = " << _xhi << std::endl ; // If list proxy does not exist, create it now _lp = new RooListProxy(Form("range::%s",GetName()),"lp",&owner,false,true) ; diff --git a/roofit/roofitcore/src/RooProdPdf.cxx b/roofit/roofitcore/src/RooProdPdf.cxx index 483510a1dd6b5..f32fd972d7d65 100644 --- a/roofit/roofitcore/src/RooProdPdf.cxx +++ b/roofit/roofitcore/src/RooProdPdf.cxx @@ -522,8 +522,6 @@ void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& int getObservablesOfCurrentPdf(pdfAllDeps, normSet); -// std::cout << GetName() << ": pdf = " << pdf->GetName() << " pdfAllDeps = " << pdfAllDeps << " pdfNSet = " << *pdfNSet << " pdfCSet = " << *pdfCSet << std::endl; - // Make list of normalization dependents for this PDF; if (!pdfNSet.empty()) { // PDF is conditional @@ -533,23 +531,18 @@ void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& int pdfNormDeps = pdfAllDeps; } -// std::cout << GetName() << ": pdfNormDeps for " << pdf->GetName() << " = " << pdfNormDeps << std::endl; - pdfIntSet.clear(); getObservablesOfCurrentPdf(pdfIntSet, intSet) ; // WVE if we have no norm deps, conditional observables should be taken out of pdfIntSet if (pdfNormDeps.empty() && !pdfCSet.empty()) { removeCommon(pdfIntSet, pdfCSet); -// std::cout << GetName() << ": have no norm deps, removing conditional observables from intset" << std::endl; } pdfIntNoNormDeps.clear(); pdfIntNoNormDeps = pdfIntSet; removeCommon(pdfIntNoNormDeps, pdfNormDeps); -// std::cout << GetName() << ": pdf = " << pdf->GetName() << " intset = " << *pdfIntSet << " pdfIntNoNormDeps = " << pdfIntNoNormDeps << std::endl; - // Check if this PDF has dependents overlapping with one of the existing terms bool done = false; int j = 0; @@ -568,8 +561,6 @@ void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& int //bool intOverlap = pdfIntSet->overlaps(*termAllDeps); if (normOverlap) { -// std::cout << GetName() << ": this term overlaps with term " << (*term) << " in normalization observables" << std::endl; - term->add(pdf); termNormDeps->add(pdfNormDeps.begin(), pdfNormDeps.end(), false); depAllList[j].add(pdfAllDeps.begin(), pdfAllDeps.end(), false); @@ -620,7 +611,6 @@ void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& int auto snap = new RooArgSet; impDeps.snapshot(*snap); factorized.imps.Add(snap); -// std::cout << GetName() << ": list of imported dependents for term " << (*term) << " set to " << impDeps << std::endl ; // Make list of cross dependents (term is self contained for these dependents, // but components import dependents from other components) @@ -628,7 +618,6 @@ void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& int snap = new RooArgSet; crossDeps->snapshot(*snap); factorized.cross.Add(snap); -// std::cout << GetName() << ": list of cross dependents for term " << (*term) << " set to " << *crossDeps << std::endl ; } } @@ -662,11 +651,6 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS const RooArgSet* iset, const char* isetRangeName) const { -// std::cout << " FOLKERT::RooProdPdf::getPartIntList(" << GetName() <<") nset = " << (nset?*nset:RooArgSet()) << std::endl -// << " _normRange = " << _normRange << std::endl -// << " iset = " << (iset?*iset:RooArgSet()) << std::endl -// << " isetRangeName = " << (isetRangeName?isetRangeName:"") << std::endl ; - // Create containers for partial integral components to be generated auto cache = std::make_unique(); @@ -684,13 +668,10 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS groupProductTerms(groupedList, outerIntDeps, factorized); // Loop over groups -// std::cout<<"FK: pdf("< ratioTerms; for (auto const& group : groupedList) { if (1 == group.size()) { -// std::cout<<"FK: Starting Single Term"<< std::endl; - RooArgSet* term = group[0]; Int_t termIdx = factorized.terms.IndexOf(term); @@ -699,25 +680,18 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS RooArgSet termNSet(*norm); RooArgSet termImpSet(*imps); - // std::cout<<"FK: termImpSet.size() = "<(term->first()), termNSet, termImpSet, normRange(), RooNameReg::str(_refRangeName)); std::ostringstream str; termImpSet.printValue(str); -// std::cout << GetName() << "inserting ratio term" << std::endl; ratioTerms[str.str()].addOwned(std::move(ratio)); } } } else { -// std::cout<<"FK: Starting Composite Term"<< std::endl; - for (auto const& term : group) { Int_t termIdx = factorized.terms.IndexOf(term); @@ -730,7 +704,6 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS // WVE we can skip this if the ref range is equal to the normalization range if (!isRangeIdentical(termNSet, _normRange, _refRangeName)) { -// std::cout << "PREPARING RATIO HERE (COMPOSITE TERM)" << std::endl ; auto ratio = makeCondPdfRatioCorr(*static_cast(term->first()), termNSet, termImpSet, normRange(), RooNameReg::str(_refRangeName)); std::ostringstream str; termImpSet.printValue(str); ratioTerms[str.str()].addOwned(std::move(ratio)); @@ -754,7 +727,6 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS // If termNset matches index of ratioTerms, insert ratio here ostringstream str; termNSet.printValue(str); if (!ratioTerms[str.str()].empty()) { -// std::cout << "MUST INSERT RATIO OBJECT IN TERM (COMPOSITE)" << *term << std::endl; term->add(ratioTerms[str.str()]); cache->_ownedList.addOwned(std::move(ratioTerms[str.str()])); } @@ -762,11 +734,7 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS } for (auto const& group : groupedList) { -// std::cout << GetName() << ":now processing group" << std::endl; -// group->Print("1"); - if (1 == group.size()) { -// std::cout << "processing atomic item" << std::endl; RooArgSet* term = group[0]; Int_t termIdx = factorized.terms.IndexOf(term); @@ -795,13 +763,11 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS cache->_denList.addOwned(std::move(func.x2)); } } else { -// std::cout << "processing composite item" << std::endl; RooArgSet compTermSet; RooArgSet compTermNorm; RooArgSet compTermNum; RooArgSet compTermDen; for (auto const &term : group) { - // std::cout << GetName() << ": processing term " << (*term) << " of composite item" << std::endl ; Int_t termIdx = factorized.terms.IndexOf(term); RooArgSet *norm = factorized.termNormDeps(termIdx); RooArgSet *integ = factorized.termIntDeps(termIdx); @@ -817,7 +783,6 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS termISet.remove(outerIntDeps, true, true); auto func = processProductTerm(nset, iset, isetRangeName, term, termNSet, termISet, true); - // std::cout << GetName() << ": created composite term component " << func.x0->GetName() << std::endl; if (func.x0) { compTermSet.add(*func.x0); if (func.isOwned) cache->_ownedList.addOwned(std::unique_ptr{func.x0}); @@ -825,15 +790,10 @@ std::unique_ptr RooProdPdf::createCacheElem(const RooArgS compTermNum.add(*func.x1.release()); compTermDen.add(*func.x2.release()); - //cache->_numList.add(*func.x1); - //cache->_denList.add(*func.x2); } } -// std::cout << GetName() << ": constructing special composite product" << std::endl; -// std::cout << GetName() << ": compTermSet = " ; compTermSet.Print("1"); - // WVE THIS NEEDS TO BE REARRANGED // compTermset is set van partial integrals to be multiplied @@ -929,8 +889,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const RooArgSet specIntDeps ; string specIntRange ; -// std::cout << "THIS IS REARRANGEPRODUCT" << std::endl ; - for (std::size_t i = 0; i < cache._partList.size(); i++) { RooAbsReal *part = static_cast(cache._partList.at(i)); @@ -938,10 +896,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const RooAbsReal *den = static_cast(cache._denList.at(i)); i++; -// std::cout << "now processing part " << part->GetName() << " of type " << part->getStringAttribute("PROD_TERM_TYPE") << std::endl ; -// std::cout << "corresponding numerator = " << num->GetName() << std::endl ; -// std::cout << "corresponding denominator = " << den->GetName() << std::endl ; - RooFormulaVar* ratio(nullptr) ; RooArgSet origNumTerm ; @@ -978,7 +932,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const func = const_cast(&static_cast(func)->integrand()); } if (func->InheritsFrom(RooProduct::Class())) { -// std::cout << "product term found: " ; func->Print() ; for(RooAbsArg * arg : static_cast(func)->components()) { if (arg->getAttribute("RATIO_TERM")) { ratio = static_cast(arg) ; @@ -989,8 +942,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const } if (ratio) { -// std::cout << "Found ratio term in numerator: " << ratio->GetName() << std::endl ; -// std::cout << "Adding only original term to numerator: " << origNumTerm << std::endl ; nomList.add(origNumTerm) ; } else { nomList.add(*num) ; @@ -1001,11 +952,9 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const for (auto iter = rangeComps.begin() ; iter != rangeComps.end() ; ++iter) { // If denominator is an integral, make a clone with the integration range adjusted to // the selected component of the normalization integral -// std::cout << "NOW PROCESSING DENOMINATOR " << den->ClassName() << "::" << den->GetName() << std::endl ; if (string("SPECINT")==part->getStringAttribute("PROD_TERM_TYPE")) { -// std::cout << "create integral: SPECINT case" << std::endl ; RooRealIntegral* orig = static_cast(num); auto specRatio = static_cast(&orig->integrand()) ; specIntDeps.add(orig->intVars()) ; @@ -1015,22 +964,11 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const //RooProduct* numtmp = (RooProduct*) specRatio->getParameter(0) ; RooProduct* dentmp = static_cast(specRatio->getParameter(1)) ; -// std::cout << "numtmp = " << numtmp->ClassName() << "::" << numtmp->GetName() << std::endl ; -// std::cout << "dentmp = " << dentmp->ClassName() << "::" << dentmp->GetName() << std::endl ; - -// std::cout << "denominator components are " << dentmp->components() << std::endl ; for (auto* parg : static_range_cast(dentmp->components())) { -// std::cout << "now processing denominator component " << parg->ClassName() << "::" << parg->GetName() << std::endl ; - if (ratio && parg->dependsOn(*ratio)) { -// std::cout << "depends in value of ratio" << std::endl ; - // Make specialize ratio instance std::unique_ptr specializedRatio{specializeRatio(*(RooFormulaVar*)ratio,iter->c_str())}; -// std::cout << "specRatio = " << std::endl ; -// specializedRatio->printComponentTree() ; - // Replace generic ratio with specialized ratio RooAbsArg *partCust(nullptr) ; if (parg->InheritsFrom(RooAddition::Class())) { @@ -1050,8 +988,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const } // Print customized denominator -// std::cout << "customized function = " << std::endl ; -// partCust->printComponentTree() ; std::unique_ptr specializedPartCust{specializeIntegral(*static_cast(partCust),iter->c_str())}; @@ -1064,14 +1000,10 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const denListList[*iter].addOwned(std::move(specIntFinal)); } else { -// std::cout << "does NOT depend on value of ratio" << std::endl ; -// parg->Print("t") ; - denListList[*iter].addOwned(specializeIntegral(*parg,iter->c_str())); } } -// std::cout << "end iteration over denominator components" << std::endl ; } else { if (ratio) { @@ -1079,7 +1011,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const std::unique_ptr specRatio{specializeRatio(*(RooFormulaVar*)ratio,iter->c_str())}; // If integral is 'Int r(y)*g(y) dy ' then divide a posteriori by r(y) -// std::cout << "have ratio, orig den = " << den->GetName() << std::endl ; RooArgSet tmp(origNumTerm) ; tmp.add(*specRatio) ; @@ -1125,9 +1056,7 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const std::unique_ptr numerator = std::make_unique(name.c_str(),name.c_str(),nomList) ; RooArgSet products ; -// std::cout << "nomList = " << nomList << std::endl ; for (map::iterator iter = denListList.begin() ; iter != denListList.end() ; ++iter) { -// std::cout << "denList[" << iter->first << "] = " << iter->second << std::endl ; name = Form("%s_denominator_comp_%s",GetName(),iter->first.c_str()) ; // WVE FIX THIS (2) RooProduct* prod_comp = new RooProduct(name.c_str(),name.c_str(),iter->second) ; @@ -1152,11 +1081,6 @@ void RooProdPdf::rearrangeProduct(RooProdPdf::CacheElem& cache) const } -// std::cout << "numerator" << std::endl ; -// numerator->printComponentTree("",0,5) ; -// std::cout << "denominator" << std::endl ; -// norm->printComponentTree("",0,5) ; - // WVE DEBUG //RooMsgService::instance().debugWorkspace()->import(RooArgSet(*numerator,*norm)) ; @@ -1193,7 +1117,6 @@ std::unique_ptr RooProdPdf::specializeIntegral(RooAbsReal& input, co // If input is integral, recreate integral but override integration range to be targetRangeName RooRealIntegral* orig = static_cast(&input) ; -// std::cout << "creating integral: integrand = " << orig->integrand().GetName() << " vars = " << orig->intVars() << " range = " << targetRangeName << std::endl ; return std::unique_ptr{orig->integrand().createIntegral(orig->intVars(),targetRangeName)}; } else if (input.InheritsFrom(RooAddition::Class())) { @@ -1201,7 +1124,6 @@ std::unique_ptr RooProdPdf::specializeIntegral(RooAbsReal& input, co // If input is sum of integrals, recreate integral from first component of set, but override integration range to be targetRangeName RooAddition* orig = static_cast(&input) ; RooRealIntegral* origInt = static_cast(orig->list1().first()) ; -// std::cout << "creating integral from addition: integrand = " << origInt->integrand().GetName() << " vars = " << origInt->intVars() << " range = " << targetRangeName << std::endl ; return std::unique_ptr{origInt->integrand().createIntegral(origInt->intVars(),targetRangeName)}; } @@ -1508,7 +1430,6 @@ double RooProdPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, co } double val = calculate(*cache,true) ; -// std::cout << "RPP::aIWN(" << GetName() << ") ,code = " << code << ", value = " << val << std::endl ; return val ; } @@ -1987,7 +1908,6 @@ void RooProdPdf::setCacheAndTrackHints(RooArgSet& trackNodes) if (parg->canNodeBeCached()==Always) { trackNodes.add(*parg) ; -// std::cout << "tracking node RooProdPdf component " << parg << " " << parg->ClassName() << "::" << parg->GetName() << std::endl ; // Additional processing to fix normalization sets in case product defines conditional observables if (RooArgSet* pdf_nset = findPdfNSet(static_cast(*parg))) { diff --git a/roofit/roofitcore/src/RooRealMPFE.cxx b/roofit/roofitcore/src/RooRealMPFE.cxx index cc3529a15adb9..33b80b33a21f0 100644 --- a/roofit/roofitcore/src/RooRealMPFE.cxx +++ b/roofit/roofitcore/src/RooRealMPFE.cxx @@ -430,13 +430,11 @@ void RooRealMPFE::calculate() const // Start asynchronous calculation of arg value if (_state==Initialize) { - // std::cout << "RooRealMPFE::calculate(" << GetName() << ") initializing" << std::endl ; const_cast(this)->initialize() ; } // Inline mode -- Calculate value now if (_state==Inline) { - // std::cout << "RooRealMPFE::calculate(" << GetName() << ") performing Inline calculation NOW" << std::endl ; _value = _arg ; clearValueDirty() ; } @@ -444,7 +442,6 @@ void RooRealMPFE::calculate() const #ifndef _WIN32 // Compare current value of variables with saved values and send changes to server if (_state==Client) { - // std::cout << "RooRealMPFE::calculate(" << GetName() << ") state is Client trigger remote calculation" << std::endl ; Int_t i(0) ; //for (i=0 ; i<_vars.size() ; i++) { @@ -468,7 +465,6 @@ void RooRealMPFE::calculate() const } if ( valChanged || constChanged || _forceCalc) { - //cout << "RooRealMPFE::calculate(" << GetName() << " variable " << var->GetName() << " changed " << std::endl ; if (_verboseClient) std::cout << "RooRealMPFE::calculate(" << GetName() << ") variable " << _vars.at(i)->GetName() << " changed" << std::endl ; if (constChanged) { @@ -535,19 +531,15 @@ double RooRealMPFE::getValV(const RooArgSet* /*nset*/) const if (isValueDirty()) { // Cache is dirty, no calculation has been started yet - //cout << "RooRealMPFE::getValF(" << GetName() << ") cache is dirty, calling calculate and evaluate" << std::endl ; calculate() ; _value = evaluate() ; } else if (_calcInProgress) { - //cout << "RooRealMPFE::getValF(" << GetName() << ") calculation in progress, calling evaluate" << std::endl ; // Cache is clean and calculation is in progress _value = evaluate() ; } else { - //cout << "RooRealMPFE::getValF(" << GetName() << ") cache is clean, doing nothing" << std::endl ; // Cache is clean and calculated value is in cache } -// std::cout << "RooRealMPFE::getValV(" << GetName() << ") value = " << Form("%5.10f",_value) << std::endl ; return _value ; } diff --git a/roofit/roofitcore/src/RooSimultaneous.cxx b/roofit/roofitcore/src/RooSimultaneous.cxx index 7169024564fde..621d3b53ecb97 100644 --- a/roofit/roofitcore/src/RooSimultaneous.cxx +++ b/roofit/roofitcore/src/RooSimultaneous.cxx @@ -55,6 +55,7 @@ in each category. #include "RooBinSamplingPdf.h" #include "RooCategory.h" #include "RooCmdConfig.h" +#include "RooCompositeDataStore.h" #include "RooDataHist.h" #include "RooDataSet.h" #include "RooGlobalFunc.h" @@ -596,6 +597,43 @@ RooPlot* RooSimultaneous::plotOn(RooPlot *frame, RooLinkedList& cmdList) const // Sanity checks if (plotSanityChecks(frame)) return frame ; + // Special case: if an asymmetry is requested with respect to our index + // category, we cannot reroute the plotting to the component pdfs. The + // component pdfs don't depend on the index category, so the asymmetry engine + // in the base class would not be able to split them by index state. Instead, + // we delegate directly to the base class implementation, which constructs the + // asymmetry from the two index-state component pdfs (see the overridden + // createAsymmetryComponent() and GitHub issue #14255). + if (auto *asymCmd = static_cast(cmdList.FindObject("Asymmetry"))) { + auto *asymCat = dynamic_cast(asymCmd->getObject(0)); + if (asymCat && asymCat == &_indexCat.arg()) { + + RooLinkedList cmdList2(cmdList); + + // The base-class asymmetry-plotting engine averages the projection over + // the projection dataset. This is not supported for the composite data + // stores that back datasets with a category index, so we flatten such a + // projection dataset into a plain (vector-backed) copy first. Both the + // copy and the replacement command must outlive the plotOn() call below, + // because the command list only stores pointers to them. + std::unique_ptr flatProjData; + RooCmdArg newProjWData; + if (auto *projWData = static_cast(cmdList2.FindObject("ProjData"))) { + auto *projData = dynamic_cast(projWData->getObject(1)); + if (projData && dynamic_cast(projData->store())) { + flatProjData = std::make_unique(projData->GetName(), projData->GetTitle(), *projData->get(), + RooFit::Import(*const_cast(projData))); + const RooArgSet *projDataSet = projWData->getSet(0); + newProjWData = projDataSet ? RooFit::ProjWData(*projDataSet, *flatProjData) + : RooFit::ProjWData(*flatProjData); + replaceOrAdd(cmdList2, newProjWData); + } + } + + return RooAbsReal::plotOn(frame, cmdList2); + } + } + // Extract projection configuration from command list RooCmdConfig pc("RooSimultaneous::plotOn(" + std::string(GetName()) + ")"); pc.defineString("sliceCatState","SliceCat",0,"",true) ; @@ -904,6 +942,32 @@ RooPlot* RooSimultaneous::plotOn(RooPlot *frame, RooLinkedList& cmdList) const } +//////////////////////////////////////////////////////////////////////////////// +/// Build the component function of an asymmetry plot (see +/// RooAbsReal::plotAsymOn()) for a fixed state of the asymmetry category. +/// +/// When the asymmetry is requested in our own index category, the component for +/// a given index state is simply the corresponding pdf. We return a clone of +/// that pdf directly instead of a RooSimultaneous with a pinned index, because +/// a RooSimultaneous compiles its per-category observables with a category +/// prefix. That prefix makes it incompatible with the vectorized evaluation +/// backend that averages the asymmetry over the projection data, and would +/// otherwise silently yield a flat (zero) asymmetry (see issue #14255). For any +/// other asymmetry category we fall back to the generic implementation. + +std::unique_ptr +RooSimultaneous::createAsymmetryComponent(const RooAbsCategoryLValue &asymCat, const RooAbsCategoryLValue &asymCatState) const +{ + if (&asymCat == &_indexCat.arg()) { + const std::string &label = _indexCat.arg().lookupName(asymCatState.getCurrentIndex()); + if (RooAbsPdf *pdf = getPdf(label)) { + return RooHelpers::cloneTreeWithSameParameters(static_cast(*pdf)); + } + } + return RooAbsReal::createAsymmetryComponent(asymCat, asymCatState); +} + + //////////////////////////////////////////////////////////////////////////////// /// Interface function used by test statistics to freeze choice of observables /// for interpretation of fraction coefficients. Needed here because a RooSimultaneous diff --git a/roofit/roofitcore/src/RooVectorDataStore.cxx b/roofit/roofitcore/src/RooVectorDataStore.cxx index 2f1b22d1b64b8..ac9edd0cdc2d5 100644 --- a/roofit/roofitcore/src/RooVectorDataStore.cxx +++ b/roofit/roofitcore/src/RooVectorDataStore.cxx @@ -800,9 +800,6 @@ void RooVectorDataStore::cacheArgs(const RooAbsArg* owner, RooArgSet& newVarSet, } // WVE need to prune tracking entries _below_ constant nodes as the're not needed -// std::cout << "Number of Cache-and-Tracked args are " << trackArgs.size() << std::endl ; -// std::cout << "Compound ordered cache parameters = " << std::endl ; -// orderedArgs.Print("v") ; checkInit() ; @@ -842,25 +839,18 @@ void RooVectorDataStore::cacheArgs(const RooAbsArg* owner, RooArgSet& newVarSet, RooArgSet* normSet(nullptr) ; const char* catNset = arg->getStringAttribute("CATNormSet") ; if (catNset) { -// std::cout << "RooVectorDataStore::cacheArgs() cached node " << arg->GetName() << " has a normalization set specification CATNormSet = " << catNset << std::endl ; - RooArgSet anset = RooHelpers::selectFromArgSet(nset ? *nset : RooArgSet{}, catNset); ownedNsets.emplace_back(anset.selectCommon(*argObs)); normSet = ownedNsets.back().get(); } const char* catCset = arg->getStringAttribute("CATCondSet") ; if (catCset) { -// std::cout << "RooVectorDataStore::cacheArgs() cached node " << arg->GetName() << " has a conditional observable set specification CATCondSet = " << catCset << std::endl ; - RooArgSet acset = RooHelpers::selectFromArgSet(nset ? *nset : RooArgSet{}, catCset); argObs->remove(acset,true,true) ; normSet = argObs ; } // now construct normalization set for component from cset/nset spec -// if (normSet) { -// std::cout << "RooVectorDaraStore::cacheArgs() component " << arg->GetName() << " has custom normalization set " << *normSet << std::endl ; -// } nsetList.push_back(normSet) ; } diff --git a/roofit/roofitcore/test/testRooSimultaneous.cxx b/roofit/roofitcore/test/testRooSimultaneous.cxx index 31948fb7d2657..56a0e674ab7b7 100644 --- a/roofit/roofitcore/test/testRooSimultaneous.cxx +++ b/roofit/roofitcore/test/testRooSimultaneous.cxx @@ -786,3 +786,45 @@ TEST(RooSimultaneous, ExpectedDataWithNonIntegerWeights) EXPECT_FLOAT_EQ(tab->get("a"), ws.var("coeff_a")->getVal()); EXPECT_FLOAT_EQ(tab->get("b"), ws.var("coeff_b")->getVal()); } + +/// GitHub issue #14255. +/// Asymmetry plots with respect to the index category of a RooSimultaneous +/// should work. The asymmetry of two Gaussians in a shared observable, sitting +/// in the +1 and -1 states of the index category, has the analytic form +/// (G+ - G-) / (G+ + G-), which we compare the plotted curve against. +TEST(RooSimultaneous, AsymmetryPlot) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING}; + + RooWorkspace ws; + ws.factory("Gaussian::gauss_A(x[-10, 10], -1.0, 1.0)"); + ws.factory("Gaussian::gauss_B(x, +1.0, 1.0)"); + ws.factory("ExtendPdf::pdf_A(gauss_A, n_A[10000.])"); + ws.factory("ExtendPdf::pdf_B(gauss_B, n_B[10000.])"); + ws.factory("SIMUL::simPdf(sample[A=-1, B=+1], A=pdf_A, B=pdf_B)"); + + RooRealVar &x = *ws.var("x"); + RooCategory &sample = *ws.cat("sample"); + + std::unique_ptr data{ws.pdf("simPdf")->generate({x, sample}, 10000)}; + + std::unique_ptr frame{x.frame()}; + // Note: the projection dataset uses a composite data store, which the plot + // must handle transparently. + ws.pdf("simPdf")->plotOn(frame.get(), Asymmetry(sample), ProjWData(sample, *data)); + + RooCurve *curve = frame->getCurve(); + ASSERT_NE(curve, nullptr); + + auto analytic = [](double xv) { + double gp = std::exp(-0.5 * (xv - 1.0) * (xv - 1.0)); + double gn = std::exp(-0.5 * (xv + 1.0) * (xv + 1.0)); + return (gp - gn) / (gp + gn); + }; + + for (double xv : {-4.0, -2.0, -1.0, 0.0, 1.0, 2.0, 4.0}) { + EXPECT_NEAR(curve->interpolate(xv), analytic(xv), 1e-6) << "at x = " << xv; + } +}