81 ClassImp(RooAbsAnaConvPdf);
87 RooAbsAnaConvPdf::RooAbsAnaConvPdf() :
99 RooAbsAnaConvPdf::RooAbsAnaConvPdf(
const char *name,
const char *title,
100 const RooResolutionModel& model, RooRealVar& cVar) :
101 RooAbsPdf(name,title), _isCopy(kFALSE),
102 _model(
"!model",
"Original resolution model",this,(RooResolutionModel&)model,kFALSE,kFALSE),
103 _convVar(
"!convVar",
"Convolution variable",this,cVar,kFALSE,kFALSE),
104 _convSet(
"!convSet",
"Set of resModel X basisFunc convolutions",this),
105 _convNormSet(nullptr),
106 _coefNormMgr(this,10),
109 _convNormSet =
new RooArgSet(cVar,
"convNormSet") ;
110 _model.absArg()->setAttribute(
"NOCacheAndTrack") ;
117 RooAbsAnaConvPdf::RooAbsAnaConvPdf(
const RooAbsAnaConvPdf& other,
const char* name) :
118 RooAbsPdf(other,name), _isCopy(kTRUE),
119 _model(
"!model",this,other._model),
120 _convVar(
"!convVar",this,other._convVar),
121 _convSet(
"!convSet",this,other._convSet),
123 _convNormSet(other._convNormSet? new RooArgSet(*other._convNormSet) : new RooArgSet() ),
124 _coefNormMgr(other._coefNormMgr,this),
125 _codeReg(other._codeReg)
128 if (_model.absArg()) {
129 _model.absArg()->setAttribute(
"NOCacheAndTrack") ;
131 other._basisList.snapshot(_basisList);
139 RooAbsAnaConvPdf::~RooAbsAnaConvPdf()
142 delete _convNormSet ;
146 std::vector<RooAbsArg*> tmp(_convSet.begin(), _convSet.end());
148 for (
auto arg : tmp) {
149 _convSet.remove(*arg) ;
169 Int_t RooAbsAnaConvPdf::declareBasis(
const char* expression,
const RooArgList& params)
173 coutE(InputArguments) <<
"RooAbsAnaConvPdf::declareBasis(" << GetName() <<
"): ERROR attempt to "
174 <<
" declare basis functions in a copied RooAbsAnaConvPdf" << endl ;
179 if (!((RooResolutionModel*)_model.absArg())->isBasisSupported(expression)) {
180 coutE(InputArguments) <<
"RooAbsAnaConvPdf::declareBasis(" << GetName() <<
"): resolution model "
181 << _model.absArg()->GetName()
182 <<
" doesn't support basis function " << expression << endl ;
187 RooArgList basisArgs(_convVar.arg()) ;
188 basisArgs.add(params) ;
190 TString basisName(expression) ;
191 TIterator* iter = basisArgs.createIterator() ;
193 while(((arg=(RooAbsArg*)iter->Next()))) {
194 basisName.Append(
"_") ;
195 basisName.Append(arg->GetName()) ;
199 RooFormulaVar* basisFunc =
new RooFormulaVar(basisName,expression,basisArgs) ;
200 basisFunc->setAttribute(
"RooWorkspace::Recycle") ;
201 basisFunc->setAttribute(
"NOCacheAndTrack") ;
202 basisFunc->setOperMode(operMode()) ;
203 _basisList.addOwned(*basisFunc) ;
206 RooAbsReal* conv = ((RooResolutionModel*)_model.absArg())->convolution(basisFunc,
this) ;
208 coutE(InputArguments) <<
"RooAbsAnaConvPdf::declareBasis(" << GetName() <<
"): unable to construct convolution with basis function '"
209 << expression <<
"'" << endl ;
212 _convSet.add(*conv) ;
214 return _convSet.index(conv) ;
222 Bool_t RooAbsAnaConvPdf::changeModel(
const RooResolutionModel& newModel)
224 RooArgList newConvSet ;
225 Bool_t allOK(kTRUE) ;
226 for (
auto convArg : _convSet) {
227 auto conv =
static_cast<RooResolutionModel*
>(convArg);
230 RooResolutionModel* newConv = newModel.convolution((RooFormulaVar*)&conv->basis(),
this) ;
231 if (!newConvSet.add(*newConv)) {
240 std::for_each(newConvSet.begin(), newConvSet.end(), [](RooAbsArg* arg){
delete arg;});
246 _convSet.removeAll() ;
247 _convSet.addOwned(newConvSet) ;
250 replaceServer((RooAbsArg&)_model.arg(),(RooAbsArg&)newModel,kFALSE,kFALSE) ;
252 _model.setArg((RooResolutionModel&)newModel) ;
266 RooAbsGenContext* RooAbsAnaConvPdf::genContext(
const RooArgSet &vars,
const RooDataSet *prototype,
267 const RooArgSet* auxProto, Bool_t verbose)
const
270 RooResolutionModel* conv =
dynamic_cast<RooResolutionModel*
>(_model.absArg());
273 RooArgSet* modelDep = _model.absArg()->getObservables(&vars) ;
274 modelDep->remove(*convVar(),kTRUE,kTRUE) ;
275 Int_t numAddDep = modelDep->getSize() ;
280 Bool_t pdfCanDir = (getGenerator(*convVar(),dummy) != 0) ;
281 Bool_t resCanDir = conv && (conv->getGenerator(*convVar(),dummy)!=0) && conv->isDirectGenSafe(*convVar()) ;
283 if (numAddDep>0 || !pdfCanDir || !resCanDir) {
287 if (numAddDep>0) reason +=
"Resolution model has more observables than the convolution variable. " ;
288 if (!pdfCanDir) reason +=
"PDF does not support internal generation of convolution observable. " ;
289 if (!resCanDir) reason +=
"Resolution model does not support internal generation of convolution observable. " ;
291 coutI(Generation) <<
"RooAbsAnaConvPdf::genContext(" << GetName() <<
") Using regular accept/reject generator for convolution p.d.f because: " << reason.c_str() << endl ;
292 return new RooGenContext(*
this,vars,prototype,auxProto,verbose) ;
295 RooAbsGenContext* context = conv->modelGenContext(*
this, vars, prototype, auxProto, verbose);
296 if (context)
return context;
299 return new RooConvGenContext(*
this,vars,prototype,auxProto,verbose) ;
309 Bool_t RooAbsAnaConvPdf::isDirectGenSafe(
const RooAbsArg& arg)
const
313 if (!TString(_convVar.absArg()->GetName()).CompareTo(arg.GetName()) &&
314 dynamic_cast<RooTruthModel*>(_model.absArg())) {
318 return RooAbsPdf::isDirectGenSafe(arg) ;
326 const RooRealVar* RooAbsAnaConvPdf::convVar()
const
328 RooResolutionModel* conv = (RooResolutionModel*) _convSet.at(0) ;
329 if (!conv)
return 0 ;
330 return &conv->convVar() ;
341 Double_t RooAbsAnaConvPdf::evaluate()
const
346 for (
auto convArg : _convSet) {
347 auto conv =
static_cast<RooAbsPdf*
>(convArg);
348 Double_t coef = coefficient(index++) ;
350 Double_t c = conv->getVal(0) ;
352 cxcoutD(Eval) <<
"RooAbsAnaConvPdf::evaluate(" << GetName() <<
") val += coef*conv [" << index-1 <<
"/"
353 << _convSet.getSize() <<
"] coef = " << r <<
" conv = " << c << endl ;
354 result += conv->getVal(0)*coef ;
356 cxcoutD(Eval) <<
"RooAbsAnaConvPdf::evaluate(" << GetName() <<
") [" << index-1 <<
"/" << _convSet.getSize() <<
"] coef = 0" << endl ;
378 Int_t RooAbsAnaConvPdf::getAnalyticalIntegralWN(RooArgSet& allVars,
379 RooArgSet& analVars,
const RooArgSet* normSet2,
const char* )
const
382 if (allVars.getSize()==0)
return 0 ;
384 if (_forceNumInt)
return 0 ;
387 RooArgSet* allDeps = getObservables(allVars) ;
388 RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;
391 RooResolutionModel *conv ;
393 RooArgSet* intSetAll =
new RooArgSet(*allDeps,
"intSetAll") ;
396 RooArgSet* intCoefSet =
new RooArgSet(
"intCoefSet") ;
397 RooArgSet* intConvSet =
new RooArgSet(
"intConvSet") ;
398 TIterator* varIter = intSetAll->createIterator() ;
399 TIterator* convIter = _convSet.createIterator() ;
401 while(((arg=(RooAbsArg*) varIter->Next()))) {
404 while(((conv=(RooResolutionModel*) convIter->Next()))) {
405 if (conv->dependsOn(*arg)) ok=kFALSE ;
409 intCoefSet->add(*arg) ;
411 intConvSet->add(*arg) ;
419 RooArgSet* normCoefSet =
new RooArgSet(
"normCoefSet") ;
420 RooArgSet* normConvSet =
new RooArgSet(
"normConvSet") ;
421 RooArgSet* normSetAll = normSet ? (
new RooArgSet(*normSet,
"normSetAll")) : 0 ;
423 varIter = normSetAll->createIterator() ;
424 while(((arg=(RooAbsArg*) varIter->Next()))) {
427 while(((conv=(RooResolutionModel*) convIter->Next()))) {
428 if (conv->dependsOn(*arg)) ok=kFALSE ;
432 normCoefSet->add(*arg) ;
434 normConvSet->add(*arg) ;
442 if (intCoefSet->getSize()==0) {
443 delete intCoefSet ; intCoefSet=0 ;
445 if (intConvSet->getSize()==0) {
446 delete intConvSet ; intConvSet=0 ;
448 if (normCoefSet->getSize()==0) {
449 delete normCoefSet ; normCoefSet=0 ;
451 if (normConvSet->getSize()==0) {
452 delete normConvSet ; normConvSet=0 ;
458 Int_t masterCode(0) ;
459 std::vector<Int_t> tmp(1, 0) ;
461 masterCode = _codeReg.store(tmp, intCoefSet, intConvSet, normCoefSet, normConvSet) + 1 ;
463 analVars.add(*allDeps) ;
465 if (normSet)
delete normSet ;
466 if (normSetAll)
delete normSetAll ;
505 Double_t RooAbsAnaConvPdf::analyticalIntegralWN(Int_t code,
const RooArgSet* normSet,
const char* rangeName)
const
510 if (code==0)
return getVal(normSet) ;
513 RooArgSet *intCoefSet, *intConvSet, *normCoefSet, *normConvSet ;
514 _codeReg.retrieve(code-1,intCoefSet,intConvSet,normCoefSet,normConvSet) ;
519 if (normCoefSet==0&&normConvSet==0) {
522 Double_t integral(0) ;
523 const TNamed *_rangeName = RooNameReg::ptr(rangeName);
524 for (
auto convArg : _convSet) {
525 auto conv =
static_cast<RooResolutionModel*
>(convArg);
526 Double_t coef = getCoefNorm(index++,intCoefSet,_rangeName) ;
529 integral += coef*(_rangeName ? conv->getNormObj(0,intConvSet,_rangeName)->getVal() : conv->getNorm(intConvSet) ) ;
530 cxcoutD(Eval) <<
"RooAbsAnaConv::aiWN(" << GetName() <<
") [" << index-1 <<
"] integral += " << conv->getNorm(intConvSet) << endl ;
539 Double_t integral(0) ;
541 const TNamed *_rangeName = RooNameReg::ptr(rangeName);
542 for (
auto convArg : _convSet) {
543 auto conv =
static_cast<RooResolutionModel*
>(convArg);
545 Double_t coefInt = getCoefNorm(index,intCoefSet,_rangeName) ;
548 Double_t term = (_rangeName ? conv->getNormObj(0,intConvSet,_rangeName)->getVal() : conv->getNorm(intConvSet) ) ;
549 integral += coefInt*term ;
552 Double_t coefNorm = getCoefNorm(index,normCoefSet) ;
555 Double_t term = conv->getNorm(normConvSet) ;
556 norm += coefNorm*term ;
561 answer = integral/norm ;
576 Int_t RooAbsAnaConvPdf::getCoefAnalyticalIntegral(Int_t , RooArgSet& , RooArgSet& ,
const char* )
const
587 Double_t RooAbsAnaConvPdf::coefAnalyticalIntegral(Int_t coef, Int_t code,
const char* )
const
589 if (code==0)
return coefficient(coef) ;
590 coutE(InputArguments) <<
"RooAbsAnaConvPdf::coefAnalyticalIntegral(" << GetName() <<
") ERROR: unrecognized integration code: " << code << endl ;
607 Bool_t RooAbsAnaConvPdf::forceAnalyticalInt(
const RooAbsArg& )
const
618 Double_t RooAbsAnaConvPdf::getCoefNorm(Int_t coefIdx,
const RooArgSet* nset,
const TNamed* rangeName)
const
620 if (nset==0)
return coefficient(coefIdx) ;
622 CacheElem* cache = (CacheElem*) _coefNormMgr.getObj(nset,0,0,rangeName) ;
625 cache =
new CacheElem ;
629 makeCoefVarList(cache->_coefVarList) ;
631 for (i=0 ; i<cache->_coefVarList.getSize() ; i++) {
632 RooAbsReal* coefInt =
static_cast<RooAbsReal&
>(*cache->_coefVarList.at(i)).createIntegral(*nset,RooNameReg::str(rangeName)) ;
633 cache->_normList.addOwned(*coefInt) ;
636 _coefNormMgr.setObj(nset,0,cache,rangeName) ;
639 return ((RooAbsReal*)cache->_normList.at(coefIdx))->getVal() ;
647 void RooAbsAnaConvPdf::makeCoefVarList(RooArgList& varList)
const
650 for (Int_t i=0 ; i<_convSet.getSize() ; i++) {
651 RooArgSet* cvars = coefVars(i) ;
652 RooAbsReal* coefVar =
new RooConvCoefVar(Form(
"%s_coefVar_%d",GetName(),i),
"coefVar",*
this,i,cvars) ;
653 varList.addOwned(*coefVar) ;
663 RooArgSet* RooAbsAnaConvPdf::coefVars(Int_t )
const
665 RooArgSet* cVars = getParameters((RooArgSet*)0) ;
666 std::vector<RooAbsArg*> tmp;
667 for (
auto arg : *cVars) {
668 for (
auto convSetArg : _convSet) {
669 if (convSetArg->dependsOn(*arg)) {
675 cVars->remove(tmp.begin(), tmp.end(),
true,
true);
689 void RooAbsAnaConvPdf::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent)
const
691 RooAbsPdf::printMultiline(os,contents,verbose,indent);
693 os << indent <<
"--- RooAbsAnaConvPdf ---" << endl;
694 TIterator* iter = _convSet.createIterator() ;
695 RooResolutionModel* conv ;
696 while (((conv=(RooResolutionModel*)iter->Next()))) {
697 conv->printMultiline(os,contents,verbose,indent) ;
704 void RooAbsAnaConvPdf::setCacheAndTrackHints(RooArgSet& trackNodes)
706 RooFIter citer = _convSet.fwdIterator() ;
708 while ((carg=citer.next())) {
709 if (carg->canNodeBeCached()==Always) {
710 trackNodes.add(*carg) ;