#include "RooFit.h"
#include "RooAcceptReject.h"
#include "RooAcceptReject.h"
#include "RooAbsReal.h"
#include "RooCategory.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooRandom.h"
#include "RooErrorHandler.h"
#include "TString.h"
#include "TIterator.h"
#include <assert.h>
ClassImp(RooAcceptReject)
;
RooAcceptReject::RooAcceptReject(const RooAbsReal &func, const RooArgSet &genVars, const RooAbsReal* maxFuncVal, Bool_t verbose) :
TNamed(func), _cloneSet(0), _funcClone(0), _funcMaxVal(maxFuncVal), _verbose(verbose)
{
RooArgSet nodes(func,func.GetName());
_cloneSet= (RooArgSet*) nodes.snapshot(kTRUE);
if (!_cloneSet) {
cout << "RooAcceptReject::RooAcceptReject(" << GetName() << ") Couldn't deep-clone function, abort," << endl ;
RooErrorHandler::softAbort() ;
}
_funcClone = (RooAbsReal*)_cloneSet->find(func.GetName());
_realSampleDim= 0;
_catSampleMult= 1;
_isValid= kTRUE;
TIterator *iterator= genVars.createIterator();
const RooAbsArg *found = 0;
const RooAbsArg *arg = 0;
while((arg= (const RooAbsArg*)iterator->Next())) {
if(arg->isDerived()) {
cout << fName << "::" << ClassName() << ": cannot generate values for derived \""
<< arg->GetName() << "\"" << endl;
_isValid= kFALSE;
continue;
}
found= (const RooAbsArg*)_cloneSet->find(arg->GetName());
if(found) {
arg= found;
const RooAbsCategory * cat = dynamic_cast<const RooAbsCategory*>(found) ;
if (cat) {
_catSampleMult *= cat->numTypes() ;
} else {
_realSampleDim++;
}
}
else {
arg= _cloneSet->addClone(*arg);
}
assert(0 != arg);
const RooCategory *catVar= dynamic_cast<const RooCategory*>(arg);
const RooRealVar *realVar= dynamic_cast<const RooRealVar*>(arg);
if(0 != catVar) {
_catVars.add(*catVar);
}
else if(0 != realVar) {
if(realVar->hasMin() && realVar->hasMax()) {
_realVars.add(*realVar);
}
else {
cout << fName << "::" << ClassName() << ": cannot generate values for \""
<< realVar->GetName() << "\" with unbound range" << endl;
_isValid= kFALSE;
}
}
else {
cout << fName << "::" << ClassName() << ": cannot generate values for \""
<< arg->GetName() << "\" with unexpected type" << endl;
_isValid= kFALSE;
}
}
delete iterator;
if(!_isValid) {
cout << fName << "::" << ClassName() << ": constructor failed with errors" << endl;
return;
}
if (!_funcMaxVal) {
if(_realSampleDim > _maxSampleDim) {
_minTrials= _minTrialsArray[_maxSampleDim]*_catSampleMult;
cout << fName << "::" << ClassName() << ": WARNING: generating " << _realSampleDim
<< " variables with accept-reject may not be accurate" << endl;
}
else {
_minTrials= _minTrialsArray[_realSampleDim]*_catSampleMult;
}
} else {
_minTrials=0 ;
}
if(_verbose) {
cout << fName << "::" << ClassName() << ":" << endl
<< " Initializing accept-reject generator for" << endl << " ";
_funcClone->Print();
if (_funcMaxVal) {
cout << " Function maximum provided, no trial sampling performed" << endl ;
} else {
cout << " Real sampling dimension is " << _realSampleDim << endl;
cout << " Category sampling multiplier is " << _catSampleMult << endl ;
cout << " Min sampling trials is " << _minTrials << endl;
}
cout << " Will generate category vars ";
TString indent(" ");
_catVars.printToStream(cout,Standard,indent);
cout << " Will generate real vars ";
_realVars.printToStream(cout,Standard,indent);
}
_funcValStore= dynamic_cast<RooRealVar*>(_funcClone->createFundamental());
assert(0 != _funcValStore);
RooArgSet cacheArgs(_catVars);
cacheArgs.add(_realVars);
cacheArgs.add(*_funcValStore);
_cache= new RooDataSet("cache","Accept-Reject Event Cache",cacheArgs);
assert(0 != _cache);
const RooArgSet *cacheVars= _cache->get();
assert(0 != cacheVars);
_funcClone->recursiveRedirectServers(*cacheVars,kFALSE);
const RooArgSet *dataVars= _cache->get();
_catVars.replace(*dataVars);
_realVars.replace(*dataVars);
_funcValPtr= (RooRealVar*)dataVars->find(_funcValStore->GetName());
_nextCatVar= _catVars.createIterator();
_nextRealVar= _realVars.createIterator();
assert(0 != _nextCatVar && 0 != _nextRealVar);
_maxFuncVal= 0;
_funcSum= 0;
_totalEvents= 0;
_eventsUsed= 0;
}
RooAcceptReject::~RooAcceptReject() {
delete _cache ;
delete _nextCatVar;
delete _nextRealVar;
delete _cloneSet;
delete _funcValStore;
}
void RooAcceptReject::printToStream(ostream &os, PrintOption , TString ) const
{
oneLinePrint(os,*this);
}
void RooAcceptReject::attachParameters(const RooArgSet& vars)
{
RooArgSet newParams(vars) ;
newParams.remove(*_cache->get(),kTRUE,kTRUE) ;
_funcClone->recursiveRedirectServers(newParams) ;
}
const RooArgSet *RooAcceptReject::generateEvent(UInt_t remaining) {
const RooArgSet *event= _cache->get();
if(event->getSize() == 1) return event;
if (!_funcMaxVal) {
while(_totalEvents < _minTrials) {
addEventToCache();
if (_cache->numEntries()>1000000) {
cout << "RooAcceptReject::generateEvent: resetting event cache" << endl ;
_cache->reset() ;
_eventsUsed = 0 ;
}
}
event= 0;
while(0 == event) {
event= nextAcceptedEvent();
if(event) break;
_cache->reset();
_eventsUsed= 0;
if(_totalEvents*_maxFuncVal <= 0) {
cout << "RooAcceptReject::generateEvent: cannot estimate efficiency...giving up" << endl;
return 0;
}
Double_t eff= _funcSum/(_totalEvents*_maxFuncVal);
Int_t extra= 1 + (Int_t)(1.05*remaining/eff);
if(_verbose) {
cout << "RooAcceptReject::generateEvent: adding " << extra << " events to the cache" << endl;
}
Double_t oldMax(_maxFuncVal);
while(extra--) addEventToCache();
if(_verbose && (_maxFuncVal > oldMax)) {
cout << "RooAcceptReject::generateEvent: estimated function maximum increased from "
<< oldMax << " to " << _maxFuncVal << endl;
}
}
if (_eventsUsed>1000000) {
_cache->reset() ;
_eventsUsed = 0 ;
}
} else {
_maxFuncVal = _funcMaxVal->getVal() ;
event = 0 ;
while(0==event) {
addEventToCache() ;
event = nextAcceptedEvent() ;
}
}
return event;
}
const RooArgSet *RooAcceptReject::nextAcceptedEvent() {
const RooArgSet *event = 0;
while((event= _cache->get(_eventsUsed))) {
_eventsUsed++ ;
Double_t r= RooRandom::uniform();
if(r*_maxFuncVal > _funcValPtr->getVal()) continue;
if(_verbose && (_eventsUsed%1000==0)) {
cerr << "RooAcceptReject: accepted event (used " << _eventsUsed << " of "
<< _cache->numEntries() << " so far)" << endl;
}
break;
}
return event;
}
void RooAcceptReject::addEventToCache() {
_nextCatVar->Reset();
RooCategory *cat = 0;
while((cat= (RooCategory*)_nextCatVar->Next())) cat->randomize();
_nextRealVar->Reset();
RooRealVar *real = 0;
while((real= (RooRealVar*)_nextRealVar->Next())) real->randomize();
Double_t val= _funcClone->getVal();
_funcValPtr->setVal(val);
if(val > _maxFuncVal) _maxFuncVal= 1.05*val;
_funcSum+= val;
_cache->fill();
_totalEvents++;
if (_verbose &&_totalEvents%10000==0) {
cerr << "RooAcceptReject: generated " << _totalEvents << " events so far." << endl ;
}
}
Double_t RooAcceptReject::getFuncMax()
{
while(_totalEvents < _minTrials) {
addEventToCache();
if (_cache->numEntries()>1000000) {
cout << "RooAcceptReject::getFuncMax: resetting event cache" << endl ;
_cache->reset() ;
_eventsUsed = 0 ;
}
}
return _maxFuncVal ;
}
const UInt_t RooAcceptReject::_minTrialsArray[] = { 100,1000,100000,10000000 };
const UInt_t RooAcceptReject::_maxSampleDim = (sizeof(_minTrialsArray) / sizeof(int)) - 1;
ROOT page - Class index - Class Hierarchy - Top of the page
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.