// Class RooProfileLL implements the profile likelihood estimator for
// a given likelihood and set of parameters of interest. The value return by
// RooProfileLL is the input likelihood nll minimized w.r.t all nuisance parameters
// (which are all parameters except for those listed in the constructor) minus
// the -log(L) of the best fit. Note that this function is slow to evaluate
// as a MIGRAD minimization step is executed for each function evaluation
// END_HTML
#include "Riostream.h"
#include "RooFit.h"
#include "RooProfileLL.h"
#include "RooAbsReal.h"
#include "RooMinuit.h"
#include "RooMinimizer.h"
#include "RooMsgService.h"
#include "RooRealVar.h"
#include "RooMsgService.h"
using namespace std ;
ClassImp(RooProfileLL)
RooProfileLL::RooProfileLL() :
RooAbsReal("RooProfileLL","RooProfileLL"),
_nll(),
_obs("paramOfInterest","Parameters of interest",this),
_par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE),
_startFromMin(kTRUE),
_minimizer(0),
_absMinValid(kFALSE),
_absMin(0),
_neval(0)
{
_piter = _par.createIterator() ;
_oiter = _obs.createIterator() ;
}
RooProfileLL::RooProfileLL(const char *name, const char *title,
RooAbsReal& nllIn, const RooArgSet& observables) :
RooAbsReal(name,title),
_nll("input","-log(L) function",this,nllIn),
_obs("paramOfInterest","Parameters of interest",this),
_par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE),
_startFromMin(kTRUE),
_minimizer(0),
_absMinValid(kFALSE),
_absMin(0),
_neval(0)
{
RooArgSet* actualObs = nllIn.getObservables(observables) ;
RooArgSet* actualPars = nllIn.getParameters(observables) ;
_obs.add(*actualObs) ;
_par.add(*actualPars) ;
delete actualObs ;
delete actualPars ;
_piter = _par.createIterator() ;
_oiter = _obs.createIterator() ;
}
RooProfileLL::RooProfileLL(const RooProfileLL& other, const char* name) :
RooAbsReal(other,name),
_nll("nll",this,other._nll),
_obs("obs",this,other._obs),
_par("par",this,other._par),
_startFromMin(other._startFromMin),
_minimizer(0),
_absMinValid(kFALSE),
_absMin(0),
_paramFixed(other._paramFixed),
_neval(0)
{
_piter = _par.createIterator() ;
_oiter = _obs.createIterator() ;
_paramAbsMin.addClone(other._paramAbsMin) ;
_obsAbsMin.addClone(other._obsAbsMin) ;
}
RooProfileLL::~RooProfileLL()
{
if (_minimizer) {
delete _minimizer ;
}
delete _piter ;
delete _oiter ;
}
const RooArgSet& RooProfileLL::bestFitParams() const
{
validateAbsMin() ;
return _paramAbsMin ;
}
const RooArgSet& RooProfileLL::bestFitObs() const
{
validateAbsMin() ;
return _obsAbsMin ;
}
RooAbsReal* RooProfileLL::createProfile(const RooArgSet& paramsOfInterest)
{
return nll().createProfile(paramsOfInterest) ;
}
void RooProfileLL::initializeMinimizer() const
{
coutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") Creating instance of MINUIT" << endl ;
Bool_t smode = RooMsgService::instance().silentMode() ;
RooMsgService::instance().setSilentMode(kTRUE) ;
_minimizer = new MINIMIZER(const_cast<RooAbsReal&>(_nll.arg())) ;
if (!smode) RooMsgService::instance().setSilentMode(kFALSE) ;
}
Double_t RooProfileLL::evaluate() const
{
if (!_minimizer) {
initializeMinimizer() ;
}
RooArgSet* obsSetOrig = (RooArgSet*) _obs.snapshot() ;
validateAbsMin() ;
const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",kTRUE) ;
ccoutP(Eval) << "." ; ccoutP(Eval).flush() ;
if (_startFromMin) {
const_cast<RooProfileLL&>(*this)._par = _paramAbsMin ;
}
_minimizer->zeroEvalCount() ;
_minimizer->migrad() ;
_neval = _minimizer->evalCounter() ;
TIterator* iter = obsSetOrig->createIterator() ;
RooRealVar* var ;
while((var=(RooRealVar*)iter->Next())) {
RooRealVar* target = (RooRealVar*) _obs.find(var->GetName()) ;
target->setVal(var->getVal()) ;
target->setConstant(var->isConstant()) ;
}
delete iter ;
delete obsSetOrig ;
return _nll - _absMin ;
}
void RooProfileLL::validateAbsMin() const
{
if (_absMinValid) {
_piter->Reset() ;
RooAbsArg* par ;
while((par=(RooAbsArg*)_piter->Next())) {
if (_paramFixed[par->GetName()] != par->isConstant()) {
cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") constant status of parameter " << par->GetName() << " has changed from "
<< (_paramFixed[par->GetName()]?"fixed":"floating") << " to " << (par->isConstant()?"fixed":"floating")
<< ", recalculating absolute minimum" << endl ;
_absMinValid = kFALSE ;
break ;
}
}
}
if (!_absMinValid) {
cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") determining minimum likelihood for current configurations w.r.t all observable" << endl ;
if (!_minimizer) {
initializeMinimizer() ;
}
RooArgSet* obsStart = (RooArgSet*) _obs.snapshot(kFALSE) ;
if (_paramAbsMin.getSize()>0) {
const_cast<RooSetProxy&>(_par).assignValueOnly(_paramAbsMin) ;
}
if (_obsAbsMin.getSize()>0) {
const_cast<RooSetProxy&>(_obs).assignValueOnly(_obsAbsMin) ;
}
const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",kFALSE) ;
_minimizer->migrad() ;
_absMin = _nll ;
_absMinValid = kTRUE ;
_paramAbsMin.removeAll() ;
RooArgSet* tmp = (RooArgSet*) _par.selectByAttrib("Constant",kFALSE) ;
_paramAbsMin.addClone(*tmp) ;
delete tmp ;
_obsAbsMin.addClone(_obs) ;
_piter->Reset() ;
RooAbsArg* par ;
while((par=(RooAbsArg*)_piter->Next())) {
_paramFixed[par->GetName()] = par->isConstant() ;
}
if (dologI(Minimization)) {
cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") minimum found at (" ;
RooAbsReal* arg ;
Bool_t first=kTRUE ;
_oiter->Reset() ;
while ((arg=(RooAbsReal*)_oiter->Next())) {
ccxcoutI(Minimization) << (first?"":", ") << arg->GetName() << "=" << arg->getVal() ;
first=kFALSE ;
}
ccxcoutI(Minimization) << ")" << endl ;
}
const_cast<RooSetProxy&>(_obs) = *obsStart ;
delete obsStart ;
}
}
Bool_t RooProfileLL::redirectServersHook(const RooAbsCollection& , Bool_t ,
Bool_t , Bool_t )
{
if (_minimizer) {
delete _minimizer ;
_minimizer = 0 ;
}
return kFALSE ;
}