#include "RooFit.h"
#include "TIterator.h"
#include "TIterator.h"
#include "TList.h"
#include "RooAddModel.h"
#include "RooRealProxy.h"
#include "RooArgList.h"
#include "RooRandom.h"
#include "RooRealConstant.h"
ClassImp(RooAddModel)
;
RooAddModel::RooAddModel(const char *name, const char *title,
const RooArgList& modelList, const RooArgList& coefList) :
RooResolutionModel(name,title,((RooResolutionModel*)modelList.at(0))->convVar()),
_nsetCache(10),
_codeReg(10),
_genReg(10),
_genThresh(0),
_isCopy(kFALSE),
_dummyProxy("dummyProxy","dummy proxy",this,(RooRealVar&)RooRealConstant::value(0)),
_modelProxyIter(_modelProxyList.MakeIterator()),
_coefProxyIter(_coefProxyList.MakeIterator())
{
if (modelList.getSize() != coefList.getSize() + 1) {
cout << "RooAddModel::ctor(" << GetName()
<< ") ERROR: number of coefficients must be one less than number of models" << endl ;
assert(0) ;
}
RooAbsArg* refConvVar = 0;
RooResolutionModel* model ;
TIterator* mIter = modelList.createIterator() ;
while((model=(RooResolutionModel*)mIter->Next())) {
if (!dynamic_cast<RooResolutionModel*>(model)) {
cout << "RooAddModel::ctor(" << GetName() << ") ERROR: " << model->GetName()
<< " is not a RooResolutionModel" << endl ;
assert(0) ;
}
if (!refConvVar) {
refConvVar = &model->convVar() ;
} else {
if (&model->convVar() != refConvVar) {
cout << "RooAddModel::ctor(" << GetName()
<< ") ERROR: models have inconsistent convolution variable" << endl ;
assert(0) ;
}
}
RooRealProxy *modelProxy = new RooRealProxy("model","model",this,*model) ;
_modelProxyList.Add(modelProxy) ;
}
delete mIter ;
RooAbsReal* coef ;
TIterator* cIter = coefList.createIterator() ;
while((coef=(RooAbsReal*)cIter->Next())) {
if (!dynamic_cast<RooAbsReal*>(coef)) {
cout << "RooAddModel::ctor(" << GetName() << ") ERROR: " << coef->GetName()
<< " is not a RooAbsReal" << endl ;
assert(0) ;
}
RooRealProxy *coefProxy = new RooRealProxy("coef","coef",this,*coef) ;
_coefProxyList.Add(coefProxy) ;
}
delete cIter ;
}
RooAddModel::RooAddModel(const RooAddModel& other, const char* name) :
RooResolutionModel(other,name),
_nsetCache(10),
_codeReg(other._codeReg),
_genReg(other._genReg),
_genThresh(0),
_isCopy(kTRUE),
_dummyProxy("dummyProxy",this,other._dummyProxy),
_modelProxyIter(_modelProxyList.MakeIterator()),
_coefProxyIter(_coefProxyList.MakeIterator())
{
TIterator *iter = other._coefProxyList.MakeIterator() ;
RooRealProxy* proxy ;
while((proxy=(RooRealProxy*)iter->Next())) {
_coefProxyList.Add(new RooRealProxy("coef",this,*proxy)) ;
}
delete iter ;
iter = other._modelProxyList.MakeIterator() ;
while((proxy=(RooRealProxy*)iter->Next())) {
_modelProxyList.Add(new RooRealProxy("model",this,*proxy)) ;
}
delete iter ;
}
RooAddModel::~RooAddModel()
{
TList ownedList ;
if (_basis && !_isCopy) {
TIterator* mIter = _modelProxyList.MakeIterator() ;
RooRealProxy* modelProxy ;
while ((modelProxy=((RooRealProxy*)mIter->Next()))) {
ownedList.Add(modelProxy->absArg()) ;
}
delete mIter ;
}
delete _coefProxyIter ;
delete _modelProxyIter ;
_coefProxyList.Delete() ;
_modelProxyList.Delete() ;
if (_basis && !_isCopy) {
ownedList.Delete() ;
}
if (_genThresh) delete[] _genThresh ;
}
RooResolutionModel* RooAddModel::convolution(RooFormulaVar* basis, RooAbsArg* owner) const
{
if (basis->findServer(0) != x.absArg()) {
cout << "RooAddModel::convolution(" << GetName()
<< ") convolution parameter of basis function and PDF don't match" << endl ;
cout << "basis->findServer(0) = " << basis->findServer(0) << " " << basis->findServer(0)->GetName() << endl ;
cout << "x.absArg() = " << x.absArg() << " " << x.absArg()->GetName() << endl ;
basis->Print("v") ;
return 0 ;
}
TString newName(GetName()) ;
newName.Append("_conv_") ;
newName.Append(basis->GetName()) ;
newName.Append("_[") ;
newName.Append(owner->GetName()) ;
newName.Append("]") ;
TString newTitle(GetTitle()) ;
newTitle.Append(" convoluted with basis function ") ;
newTitle.Append(basis->GetName()) ;
_modelProxyIter->Reset() ;
RooRealProxy* model ;
RooArgList modelList ;
while((model = (RooRealProxy*)_modelProxyIter->Next())) {
RooResolutionModel* conv = ((RooResolutionModel*)(model->absArg()))->convolution(basis,owner) ;
modelList.add(*conv) ;
}
_coefProxyIter->Reset() ;
RooRealProxy* coef ;
RooArgList coefList ;
while((coef = (RooRealProxy*)_coefProxyIter->Next())) {
coefList.add(coef->arg()) ;
}
RooAddModel* convSum = new RooAddModel(newName,newTitle,modelList,coefList) ;
convSum->changeBasis(basis) ;
return convSum ;
}
Int_t RooAddModel::basisCode(const char* name) const
{
TIterator* mIter = _modelProxyList.MakeIterator() ;
RooRealProxy* model ;
Bool_t first(kTRUE), code(0) ;
while((model = (RooRealProxy*)mIter->Next())) {
Int_t subCode = ((RooResolutionModel&)model->arg()).basisCode(name) ;
if (first) {
code = subCode ;
} else if (subCode==0) {
code = 0 ;
}
}
delete mIter ;
return code ;
}
Double_t RooAddModel::evaluate() const
{
_coefProxyIter->Reset() ;
_modelProxyIter->Reset() ;
Double_t value(0) ;
Double_t lastCoef(1) ;
const RooArgSet* nset = _dummyProxy.nset() ;
RooRealProxy* coef ;
RooRealProxy* model ;
while((coef=(RooRealProxy*)_coefProxyIter->Next())) {
model = (RooRealProxy*)_modelProxyIter->Next() ;
Double_t coefVal = coef->arg().getVal(nset) ;
if (coefVal) {
if (((RooAbsPdf&)model->arg()).isSelectedComp()) {
value += model->arg().getVal(nset)*coef->arg().getVal(nset) ;
}
lastCoef -= coef->arg().getVal(nset) ;
}
}
model = (RooRealProxy*) _modelProxyIter->Next() ;
if (((RooAbsPdf&)model->arg()).isSelectedComp()) {
value += model->arg().getVal(nset)*lastCoef ;
}
if ((lastCoef<0.0 || lastCoef>1.0) && ++_errorCount<=10) {
cout << "RooAddModel::evaluate(" << GetName()
<< " WARNING: sum of model coefficients not in range [0-1], value="
<< 1-lastCoef << endl ;
if(_errorCount == 10) cout << "(no more will be printed) ";
}
return value ;
}
Double_t RooAddModel::getNorm(const RooArgSet* nset) const
{
if (!_basis) return RooAbsPdf::getNorm(nset) ;
_coefProxyIter->Reset() ;
_modelProxyIter->Reset() ;
Double_t norm(0) ;
Double_t lastCoef(1) ;
RooRealProxy* coef ;
RooResolutionModel* model ;
while((coef=(RooRealProxy*)_coefProxyIter->Next())) {
model = (RooResolutionModel*)((RooRealProxy*)_modelProxyIter->Next())->absArg() ;
if (_verboseEval>1) {
cout << "RooAddModel::getNorm(" << GetName() << "): norm x coef = "
<< model->getNorm(nset) << " x " << (*coef) << " = "
<< model->getNorm(nset)*(*coef) ;
cout << "nset = " ; nset->Print("1") ;
}
Double_t coefVal = *coef ;
if (coefVal) {
norm += model->getNorm(nset)*(*coef) ;
lastCoef -= (*coef) ;
}
}
model = (RooResolutionModel*)((RooRealProxy*)_modelProxyIter->Next())->absArg() ;
norm += model->getNorm(nset)*lastCoef ;
if (_verboseEval>1) cout << "RooAddModel::getNorm(" << GetName() << "): norm x coef = "
<< model->getNorm(nset) << " x " << lastCoef << " = "
<< model->getNorm(nset)*lastCoef << endl ;
if ((lastCoef<0 || lastCoef>1) && ++_errorCount<=10) {
cout << "RooAddModel::evaluate(" << GetName()
<< " WARNING: sum of model coefficients not in range [0-1], value="
<< 1-lastCoef << endl ;
}
return norm ;
}
Bool_t RooAddModel::checkObservables(const RooArgSet* set) const
{
Bool_t ret(kFALSE) ;
TIterator *pIter = _modelProxyList.MakeIterator() ;
TIterator *cIter = _coefProxyList.MakeIterator() ;
RooRealProxy* coef ;
RooRealProxy* model ;
while((coef=(RooRealProxy*)cIter->Next())) {
model = (RooRealProxy*)pIter->Next() ;
if (model->arg().observableOverlaps(set,coef->arg())) {
cout << "RooAddModel::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->arg().GetName()
<< " and model " << model->arg().GetName() << " have one or more dependents in common" << endl ;
ret = kTRUE ;
}
}
delete pIter ;
delete cIter ;
return ret ;
}
void RooAddModel::normLeafServerList(RooArgSet& list) const
{
TIterator *pIter = _modelProxyList.MakeIterator() ;
RooRealProxy* proxy ;
RooResolutionModel* model ;
while((proxy = (RooRealProxy*) pIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
if (model->_norm==0) {
model->syncNormalization(proxy->nset()) ;
}
model->_norm->leafNodeServerList(&list) ;
}
delete pIter ;
}
Bool_t RooAddModel::syncNormalization(const RooArgSet* nset, Bool_t adjustProxies) const
{
if (!_nsetCache.autoCache(this,nset)) return kFALSE ;
if (_verboseEval>0) cout << "RooAddModel:syncNormalization(" << GetName()
<< ") forwarding sync request to components" << endl ;
((RooAbsPdf*) this)->setProxyNormSet(nset) ;
TIterator *pIter = _modelProxyList.MakeIterator() ;
RooRealProxy* proxy ;
RooResolutionModel* model ;
while((proxy = (RooRealProxy*)pIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
model->syncNormalization(nset,adjustProxies) ;
}
delete pIter ;
if (_basisCode==0) {
if (_verboseEval>0) {
cout << "RooAddModel::syncNormalization(" << GetName()
<< ") creating unit normalization object" << endl ;
}
TString nname(GetName()) ; nname.Append("Norm") ;
TString ntitle(GetTitle()) ; ntitle.Append(" Unit Normalization") ;
if (_norm) delete _norm ;
_norm = new RooRealVar(nname.Data(),ntitle.Data(),1) ;
}
return kTRUE ;
}
Bool_t RooAddModel::forceAnalyticalInt(const RooAbsArg& ) const
{
return (_basisCode==0) ;
}
Int_t RooAddModel::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars,
const RooArgSet* normSet, const char* ) const
{
if (_basisCode!=0) {
return 0 ;
}
_modelProxyIter->Reset() ;
RooResolutionModel* model ;
RooArgSet allAnalVars(allVars) ;
TIterator* avIter = allVars.createIterator() ;
Int_t n(0) ;
RooRealProxy* proxy ;
while((proxy=(RooRealProxy*)_modelProxyIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
RooArgSet subAnalVars ;
model->getAnalyticalIntegralWN(allVars,subAnalVars,normSet) ;
avIter->Reset() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)avIter->Next())) {
if (!subAnalVars.find(arg->GetName())) {
allAnalVars.remove(*arg,kTRUE) ;
}
}
n++ ;
}
if (allAnalVars.getSize()==0) {
delete avIter ;
return 0 ;
}
_modelProxyIter->Reset() ;
n=0 ;
Int_t* subCode = new Int_t[_modelProxyList.GetSize()] ;
Bool_t allOK(kTRUE) ;
while((proxy=(RooRealProxy*)_modelProxyIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
RooArgSet subAnalVars ;
subCode[n] = model->getAnalyticalIntegralWN(allAnalVars,subAnalVars,normSet) ;
if (subCode[n]==0) {
cout << "RooAddModel::getAnalyticalIntegral(" << GetName() << ") WARNING: component model " << model->GetName()
<< " advertises inconsistent set of integrals (e.g. (X,Y) but not X or Y individually."
<< " Distributed analytical integration disabled. Please fix model" << endl ;
allOK = kFALSE ;
}
n++ ;
}
if (!allOK) return 0 ;
analVars.add(allAnalVars) ;
Int_t masterCode = _codeReg.store(subCode,_modelProxyList.GetSize())+1 ;
delete[] subCode ;
delete avIter ;
return masterCode ;
}
Double_t RooAddModel::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* ) const
{
if (code==0) return getVal() ;
const Int_t* subCode = _codeReg.retrieve(code-1) ;
if (!subCode) {
cout << "RooAddModel::analyticalIntegral(" << GetName() << "): ERROR unrecognized integration code, " << code << endl ;
assert(0) ;
}
Double_t value(0) ;
_modelProxyIter->Reset() ;
_coefProxyIter->Reset() ;
RooAbsReal* coef ;
RooResolutionModel* model ;
Int_t i(0) ;
Double_t lastCoef(1) ;
RooRealProxy* proxy ;
while((proxy=(RooRealProxy*)_coefProxyIter->Next())) {
coef = (RooAbsReal*) proxy->absArg() ;
model = (RooResolutionModel*)((RooRealProxy*)_modelProxyIter->Next())->absArg() ;
Double_t coefVal = coef->getVal(normSet) ;
value += model->analyticalIntegralWN(subCode[i],normSet)*coefVal ;
lastCoef -= coefVal ;
i++ ;
}
model = (RooResolutionModel*) ((RooRealProxy*)_modelProxyIter->Next())->absArg() ;
value += model->analyticalIntegralWN(subCode[i],normSet)*lastCoef ;
if ((lastCoef<0 || lastCoef>1) && ++_errorCount<=10) {
cout << "RooAddModel::analyticalIntegral(" << GetName()
<< " WARNING: sum of model coefficients not in range [0-1], value="
<< 1-lastCoef << endl ;
}
return value ;
}
Bool_t RooAddModel::isDirectGenSafe(const RooAbsArg& arg) const
{
RooRealProxy* proxy ;
RooResolutionModel* model ;
_modelProxyIter->Reset() ;
while((proxy=(RooRealProxy*)_modelProxyIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
if (!model->isDirectGenSafe(arg)) return kFALSE ;
}
return kTRUE ;
}
Int_t RooAddModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
{
Int_t* subcode = new Int_t[_modelProxyList.GetSize()] ;
RooRealProxy* proxy ;
Int_t n(0) ;
RooResolutionModel* model ;
_modelProxyIter->Reset() ;
while((proxy=(RooRealProxy*)_modelProxyIter->Next())) {
model = (RooResolutionModel*) proxy->absArg() ;
RooArgSet subGenVars ;
subcode[n] = model->getGenerator(directVars,subGenVars,staticInitOK) ;
if (subcode[n]==0 || (subGenVars.getSize() != directVars.getSize())) {
return 0 ;
}
n++ ;
}
generateVars.add(directVars) ;
Int_t masterCode = _genReg.store(subcode,_modelProxyList.GetSize())+1 ;
return masterCode ;
}
void RooAddModel::initGenerator(Int_t code)
{
if (_genThresh) delete[] _genThresh ;
_genThresh = new Double_t[_modelProxyList.GetSize()+1] ;
Int_t nComp = _modelProxyList.GetSize() ;
_genThresh = new Double_t[nComp+1] ;
Int_t i=1 ;
_genThresh[0] = 0 ;
_modelProxyIter->Reset() ;
RooRealProxy* proxy ;
while((proxy=(RooRealProxy*)_modelProxyIter->Next())) {
((RooResolutionModel*)proxy->absArg())->initGenerator(code) ;
RooRealProxy *coefProxy = (RooRealProxy*) _coefProxyList.At(i-1) ;
RooAbsReal* coef = (RooAbsReal*) (coefProxy ? coefProxy->absArg() : 0) ;
if (coef) {
_genThresh[i] = coef->getVal() ;
if (i>0) _genThresh[i] += _genThresh[i-1] ;
} else {
_genThresh[i] = 1.0 ;
}
i++ ;
}
_genSubCode = _genReg.retrieve(code-1) ;
}
void RooAddModel::generateEvent(Int_t )
{
Double_t rand = RooRandom::uniform() ;
Int_t i=0 ;
Int_t nComp = _modelProxyList.GetSize() ;
for (i=0 ; i<nComp ; i++) {
if (rand>_genThresh[i] && rand<_genThresh[i+1]) {
RooResolutionModel* model = (RooResolutionModel*) ((RooRealProxy*)_modelProxyList.At(i))->absArg() ;
model->generateEvent(_genSubCode[i]) ;
return ;
}
}
}
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.