00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 #include "RooFit.h"
00059 #include "RooMsgService.h"
00060
00061 #include "Riostream.h"
00062 #include "Riostream.h"
00063 #include "RooAbsAnaConvPdf.h"
00064 #include "RooResolutionModel.h"
00065 #include "RooRealVar.h"
00066 #include "RooFormulaVar.h"
00067 #include "RooConvGenContext.h"
00068 #include "RooGenContext.h"
00069 #include "RooTruthModel.h"
00070 #include "RooConvCoefVar.h"
00071 #include "RooNameReg.h"
00072
00073 ClassImp(RooAbsAnaConvPdf)
00074 ;
00075
00076
00077
00078 RooAbsAnaConvPdf::RooAbsAnaConvPdf() :
00079 _isCopy(kFALSE),
00080 _model(0),
00081 _convVar(0),
00082 _convNormSet(0),
00083 _convSetIter(_convSet.createIterator())
00084 {
00085
00086 }
00087
00088
00089
00090
00091 RooAbsAnaConvPdf::RooAbsAnaConvPdf(const char *name, const char *title,
00092 const RooResolutionModel& model, RooRealVar& cVar) :
00093 RooAbsPdf(name,title), _isCopy(kFALSE),
00094 _model((RooResolutionModel*)&model), _convVar((RooRealVar*)&cVar),
00095 _convSet("!convSet","Set of resModel X basisFunc convolutions",this),
00096 _convNormSet(0), _convSetIter(_convSet.createIterator()),
00097 _coefNormMgr(this,10),
00098 _codeReg(10)
00099 {
00100
00101
00102
00103 _convNormSet = new RooArgSet(cVar,"convNormSet") ;
00104 }
00105
00106
00107
00108
00109 RooAbsAnaConvPdf::RooAbsAnaConvPdf(const RooAbsAnaConvPdf& other, const char* name) :
00110 RooAbsPdf(other,name), _isCopy(kTRUE),
00111 _model(other._model),
00112 _convVar(other._convVar),
00113 _convSet("!convSet",this,other._convSet),
00114 _basisList(other._basisList),
00115 _convNormSet(other._convNormSet? new RooArgSet(*other._convNormSet) : new RooArgSet() ),
00116 _convSetIter(_convSet.createIterator()),
00117 _coefNormMgr(other._coefNormMgr,this),
00118 _codeReg(other._codeReg)
00119 {
00120
00121
00122 }
00123
00124
00125
00126
00127 RooAbsAnaConvPdf::~RooAbsAnaConvPdf()
00128 {
00129
00130
00131 if (_convNormSet) {
00132 delete _convNormSet ;
00133 }
00134
00135 delete _convSetIter ;
00136
00137 if (!_isCopy) {
00138 TIterator* iter = _convSet.createIterator() ;
00139 RooAbsArg* arg ;
00140 while (((arg = (RooAbsArg*)iter->Next()))) {
00141 _convSet.remove(*arg) ;
00142 delete arg ;
00143 }
00144 delete iter ;
00145 }
00146
00147 }
00148
00149
00150
00151 Int_t RooAbsAnaConvPdf::declareBasis(const char* expression, const RooArgList& params)
00152 {
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 if (_isCopy) {
00166 coutE(InputArguments) << "RooAbsAnaConvPdf::declareBasis(" << GetName() << "): ERROR attempt to "
00167 << " declare basis functions in a copied RooAbsAnaConvPdf" << endl ;
00168 return -1 ;
00169 }
00170
00171
00172 if (!_model->isBasisSupported(expression)) {
00173 coutE(InputArguments) << "RooAbsAnaConvPdf::declareBasis(" << GetName() << "): resolution model "
00174 << _model->GetName()
00175 << " doesn't support basis function " << expression << endl ;
00176 return -1 ;
00177 }
00178
00179
00180 RooArgList basisArgs(*_convVar) ;
00181 basisArgs.add(params) ;
00182
00183 TString basisName(expression) ;
00184 TIterator* iter = basisArgs.createIterator() ;
00185 RooAbsArg* arg ;
00186 while(((arg=(RooAbsArg*)iter->Next()))) {
00187 basisName.Append("_") ;
00188 basisName.Append(arg->GetName()) ;
00189 }
00190 delete iter ;
00191
00192 RooFormulaVar* basisFunc = new RooFormulaVar(basisName,expression,basisArgs) ;
00193 basisFunc->setAttribute("RooWorkspace::Recycle") ;
00194 basisFunc->setOperMode(operMode()) ;
00195 _basisList.addOwned(*basisFunc) ;
00196
00197
00198 RooAbsReal* conv = _model->convolution(basisFunc,this) ;
00199 if (!conv) {
00200 coutE(InputArguments) << "RooAbsAnaConvPdf::declareBasis(" << GetName() << "): unable to construct convolution with basis function '"
00201 << expression << "'" << endl ;
00202 return -1 ;
00203 }
00204 _convSet.add(*conv) ;
00205
00206 return _convSet.index(conv) ;
00207 }
00208
00209
00210
00211
00212 Bool_t RooAbsAnaConvPdf::changeModel(const RooResolutionModel& newModel)
00213 {
00214
00215
00216 TIterator* cIter = _convSet.createIterator() ;
00217 RooResolutionModel* conv ;
00218 RooArgList newConvSet ;
00219 Bool_t allOK(kTRUE) ;
00220 while(((conv=(RooResolutionModel*)cIter->Next()))) {
00221
00222
00223 RooResolutionModel* newConv = newModel.convolution((RooFormulaVar*)&conv->basis(),this) ;
00224 if (!newConvSet.add(*newConv)) {
00225 allOK = kFALSE ;
00226 break ;
00227 }
00228 }
00229 delete cIter ;
00230
00231
00232 if (!allOK) {
00233
00234 TIterator* iter = newConvSet.createIterator() ;
00235 while(((conv=(RooResolutionModel*)iter->Next()))) delete conv ;
00236 delete iter ;
00237
00238 return kTRUE ;
00239 }
00240
00241
00242 _convSet.removeAll() ;
00243 _convSet.addOwned(newConvSet) ;
00244
00245 _model = (RooResolutionModel*) &newModel ;
00246 return kFALSE ;
00247 }
00248
00249
00250
00251
00252
00253 RooAbsGenContext* RooAbsAnaConvPdf::genContext(const RooArgSet &vars, const RooDataSet *prototype,
00254 const RooArgSet* auxProto, Bool_t verbose) const
00255 {
00256
00257
00258
00259
00260
00261
00262
00263 RooArgSet* modelDep = _model->getObservables(&vars) ;
00264 modelDep->remove(*convVar(),kTRUE,kTRUE) ;
00265 Int_t numAddDep = modelDep->getSize() ;
00266 delete modelDep ;
00267
00268 if (dynamic_cast<RooTruthModel*>(_model)) {
00269
00270 RooArgSet forceDirect(*convVar()) ;
00271 return new RooGenContext(*this,vars,prototype,auxProto,verbose,&forceDirect) ;
00272 }
00273
00274
00275 RooArgSet dummy ;
00276 Bool_t pdfCanDir = (getGenerator(*convVar(),dummy) != 0) ;
00277 RooResolutionModel* conv = (RooResolutionModel*) _convSet.at(0) ;
00278 Bool_t resCanDir = conv && (conv->getGenerator(*convVar(),dummy)!=0) && conv->isDirectGenSafe(*convVar()) ;
00279
00280 if (numAddDep>0 || !pdfCanDir || !resCanDir) {
00281
00282
00283 string reason ;
00284 if (numAddDep>0) reason += "Resolution model has more onservables that the convolution variable. " ;
00285 if (!pdfCanDir) reason += "PDF does not support internal generation of convolution observable. " ;
00286 if (!resCanDir) reason += "Resolution model does not support internal generation of convolution observable. " ;
00287
00288 coutI(Generation) << "RooAbsAnaConvPdf::genContext(" << GetName() << ") Using regular accept/reject generator for convolution p.d.f because: " << reason.c_str() << endl ;
00289 return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
00290 }
00291
00292
00293 return new RooConvGenContext(*this,vars,prototype,auxProto,verbose) ;
00294 }
00295
00296
00297
00298
00299 Bool_t RooAbsAnaConvPdf::isDirectGenSafe(const RooAbsArg& arg) const
00300 {
00301
00302
00303
00304
00305
00306
00307 if (!TString(_convVar->GetName()).CompareTo(arg.GetName()) &&
00308 dynamic_cast<RooTruthModel*>(_model)) {
00309 return kTRUE ;
00310 }
00311
00312 return RooAbsPdf::isDirectGenSafe(arg) ;
00313 }
00314
00315
00316
00317
00318 const RooRealVar* RooAbsAnaConvPdf::convVar() const
00319 {
00320
00321
00322 RooResolutionModel* conv = (RooResolutionModel*) _convSet.at(0) ;
00323 if (!conv) return 0 ;
00324 return &conv->convVar() ;
00325 }
00326
00327
00328
00329
00330 Double_t RooAbsAnaConvPdf::evaluate() const
00331 {
00332
00333
00334
00335
00336 Double_t result(0) ;
00337
00338 _convSetIter->Reset() ;
00339 RooAbsPdf* conv ;
00340 Int_t index(0) ;
00341 while(((conv=(RooAbsPdf*)_convSetIter->Next()))) {
00342 Double_t coef = coefficient(index++) ;
00343 if (coef!=0.) {
00344 Double_t c = conv->getVal(0) ;
00345 Double_t r = coef ;
00346 cxcoutD(Eval) << "RooAbsAnaConvPdf::evaluate(" << GetName() << ") val += coef*conv [" << index-1 << "/"
00347 << _convSet.getSize() << "] coef = " << r << " conv = " << c << endl ;
00348 result += conv->getVal(0)*coef ;
00349 } else {
00350 cxcoutD(Eval) << "RooAbsAnaConvPdf::evaluate(" << GetName() << ") [" << index-1 << "/" << _convSet.getSize() << "] coef = 0" << endl ;
00351 }
00352 }
00353
00354 return result ;
00355 }
00356
00357
00358
00359
00360 Int_t RooAbsAnaConvPdf::getAnalyticalIntegralWN(RooArgSet& allVars,
00361 RooArgSet& analVars, const RooArgSet* normSet2, const char* ) const
00362 {
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 if (allVars.getSize()==0) return 0 ;
00377
00378 if (_forceNumInt) return 0 ;
00379
00380
00381 RooArgSet* allDeps = getObservables(allVars) ;
00382 RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;
00383
00384 RooAbsArg *arg ;
00385 RooResolutionModel *conv ;
00386
00387 RooArgSet* intSetAll = new RooArgSet(*allDeps,"intSetAll") ;
00388
00389
00390 RooArgSet* intCoefSet = new RooArgSet("intCoefSet") ;
00391 RooArgSet* intConvSet = new RooArgSet("intConvSet") ;
00392 TIterator* varIter = intSetAll->createIterator() ;
00393 TIterator* convIter = _convSet.createIterator() ;
00394
00395 while(((arg=(RooAbsArg*) varIter->Next()))) {
00396 Bool_t ok(kTRUE) ;
00397 convIter->Reset() ;
00398 while(((conv=(RooResolutionModel*) convIter->Next()))) {
00399 if (conv->dependsOn(*arg)) ok=kFALSE ;
00400 }
00401
00402 if (ok) {
00403 intCoefSet->add(*arg) ;
00404 } else {
00405 intConvSet->add(*arg) ;
00406 }
00407
00408 }
00409 delete varIter ;
00410
00411
00412
00413 RooArgSet* normCoefSet = new RooArgSet("normCoefSet") ;
00414 RooArgSet* normConvSet = new RooArgSet("normConvSet") ;
00415 RooArgSet* normSetAll = normSet ? (new RooArgSet(*normSet,"normSetAll")) : 0 ;
00416 if (normSetAll) {
00417 varIter = normSetAll->createIterator() ;
00418 while(((arg=(RooAbsArg*) varIter->Next()))) {
00419 Bool_t ok(kTRUE) ;
00420 convIter->Reset() ;
00421 while(((conv=(RooResolutionModel*) convIter->Next()))) {
00422 if (conv->dependsOn(*arg)) ok=kFALSE ;
00423 }
00424
00425 if (ok) {
00426 normCoefSet->add(*arg) ;
00427 } else {
00428 normConvSet->add(*arg) ;
00429 }
00430
00431 }
00432 delete varIter ;
00433 }
00434 delete convIter ;
00435
00436 if (intCoefSet->getSize()==0) {
00437 delete intCoefSet ; intCoefSet=0 ;
00438 }
00439 if (intConvSet->getSize()==0) {
00440 delete intConvSet ; intConvSet=0 ;
00441 }
00442 if (normCoefSet->getSize()==0) {
00443 delete normCoefSet ; normCoefSet=0 ;
00444 }
00445 if (normConvSet->getSize()==0) {
00446 delete normConvSet ; normConvSet=0 ;
00447 }
00448
00449
00450
00451 Int_t masterCode(0) ;
00452 Int_t tmp(0) ;
00453
00454 masterCode = _codeReg.store(&tmp,1,intCoefSet,intConvSet,normCoefSet,normConvSet)+1 ;
00455
00456 analVars.add(*allDeps) ;
00457 delete allDeps ;
00458 if (normSet) delete normSet ;
00459 if (normSetAll) delete normSetAll ;
00460 delete intSetAll ;
00461
00462
00463
00464 return masterCode ;
00465 }
00466
00467
00468
00469
00470
00471 Double_t RooAbsAnaConvPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
00472 {
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500 if (code==0) return getVal(normSet) ;
00501
00502
00503 RooArgSet *intCoefSet, *intConvSet, *normCoefSet, *normConvSet ;
00504 _codeReg.retrieve(code-1,intCoefSet,intConvSet,normCoefSet,normConvSet) ;
00505
00506 RooResolutionModel* conv ;
00507 Int_t index(0) ;
00508 Double_t answer(0) ;
00509 _convSetIter->Reset() ;
00510
00511 if (normCoefSet==0&&normConvSet==0) {
00512
00513
00514 Double_t integral(0) ;
00515 while(((conv=(RooResolutionModel*)_convSetIter->Next()))) {
00516 Double_t coef = getCoefNorm(index++,intCoefSet,rangeName) ;
00517
00518 if (coef!=0) {
00519 integral += coef*(rangeName ? conv->getNormObj(0,intConvSet,RooNameReg::ptr(rangeName))->getVal() : conv->getNorm(intConvSet) ) ;
00520 cxcoutD(Eval) << "RooAbsAnaConv::aiWN(" << GetName() << ") [" << index-1 << "] integral += " << conv->getNorm(intConvSet) << endl ;
00521 }
00522
00523 }
00524 answer = integral ;
00525
00526 } else {
00527
00528
00529 Double_t integral(0) ;
00530 Double_t norm(0) ;
00531 while(((conv=(RooResolutionModel*)_convSetIter->Next()))) {
00532
00533 Double_t coefInt = getCoefNorm(index,intCoefSet,rangeName) ;
00534 Double_t term = (rangeName ? conv->getNormObj(0,intConvSet,RooNameReg::ptr(rangeName))->getVal() : conv->getNorm(intConvSet) ) ;
00535
00536 if (coefInt!=0) integral += coefInt*term ;
00537
00538 Double_t coefNorm = getCoefNorm(index,normCoefSet) ;
00539 term = conv->getNorm(normConvSet) ;
00540
00541 if (coefNorm!=0) norm += coefNorm*term ;
00542
00543 index++ ;
00544 }
00545 answer = integral/norm ;
00546 }
00547
00548 return answer ;
00549 }
00550
00551
00552
00553
00554 Int_t RooAbsAnaConvPdf::getCoefAnalyticalIntegral(Int_t , RooArgSet& , RooArgSet& , const char* ) const
00555 {
00556
00557
00558
00559
00560
00561
00562 return 0 ;
00563 }
00564
00565
00566
00567
00568 Double_t RooAbsAnaConvPdf::coefAnalyticalIntegral(Int_t coef, Int_t code, const char* ) const
00569 {
00570
00571
00572
00573 if (code==0) return coefficient(coef) ;
00574 coutE(InputArguments) << "RooAbsAnaConvPdf::coefAnalyticalIntegral(" << GetName() << ") ERROR: unrecognized integration code: " << code << endl ;
00575 assert(0) ;
00576 return 1 ;
00577 }
00578
00579
00580
00581
00582 Bool_t RooAbsAnaConvPdf::forceAnalyticalInt(const RooAbsArg& ) const
00583 {
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593 return kTRUE ;
00594 }
00595
00596
00597
00598
00599 Double_t RooAbsAnaConvPdf::getCoefNorm(Int_t coefIdx, const RooArgSet* nset, const char* rangeName) const
00600 {
00601
00602
00603
00604 if (nset==0) return coefficient(coefIdx) ;
00605
00606 CacheElem* cache = (CacheElem*) _coefNormMgr.getObj(nset,0,0,RooNameReg::ptr(rangeName)) ;
00607 if (!cache) {
00608
00609 cache = new CacheElem ;
00610
00611
00612 Int_t i ;
00613 makeCoefVarList(cache->_coefVarList) ;
00614
00615 for (i=0 ; i<cache->_coefVarList.getSize() ; i++) {
00616 RooAbsReal* coefInt = static_cast<RooAbsReal&>(*cache->_coefVarList.at(i)).createIntegral(*nset,rangeName) ;
00617 cache->_normList.addOwned(*coefInt) ;
00618 }
00619
00620 _coefNormMgr.setObj(nset,0,cache,RooNameReg::ptr(rangeName)) ;
00621 }
00622
00623 return ((RooAbsReal*)cache->_normList.at(coefIdx))->getVal() ;
00624 }
00625
00626
00627
00628
00629 void RooAbsAnaConvPdf::makeCoefVarList(RooArgList& varList) const
00630 {
00631
00632
00633
00634 for (Int_t i=0 ; i<_convSet.getSize() ; i++) {
00635 RooArgSet* cvars = coefVars(i) ;
00636 RooAbsReal* coefVar = new RooConvCoefVar(Form("%s_coefVar_%d",GetName(),i),"coefVar",*this,i,cvars) ;
00637 varList.addOwned(*coefVar) ;
00638 delete cvars ;
00639 }
00640
00641 }
00642
00643
00644
00645 RooArgSet* RooAbsAnaConvPdf::coefVars(Int_t ) const
00646 {
00647
00648
00649 RooArgSet* cVars = getParameters((RooArgSet*)0) ;
00650 TIterator* iter = cVars->createIterator() ;
00651 RooAbsArg* arg ;
00652 Int_t i ;
00653 while(((arg=(RooAbsArg*)iter->Next()))) {
00654 for (i=0 ; i<_convSet.getSize() ; i++) {
00655 if (_convSet.at(i)->dependsOn(*arg)) {
00656 cVars->remove(*arg,kTRUE) ;
00657 }
00658 }
00659 }
00660 delete iter ;
00661 return cVars ;
00662 }
00663
00664
00665
00666
00667
00668 void RooAbsAnaConvPdf::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
00669 {
00670
00671
00672
00673
00674
00675 RooAbsPdf::printMultiline(os,contents,verbose,indent);
00676
00677 os << indent << "--- RooAbsAnaConvPdf ---" << endl;
00678 TIterator* iter = _convSet.createIterator() ;
00679 RooResolutionModel* conv ;
00680 while (((conv=(RooResolutionModel*)iter->Next()))) {
00681 conv->printMultiline(os,contents,verbose,indent) ;
00682 }
00683 }
00684
00685