// RooMinimizer is a wrapper class around ROOT::Fit:Fitter that
// provides a seamless interface between the minimizer functionality
// and the native RooFit interface.
// <p>
// By default the Minimizer is MINUIT.
// <p>
// RooMinimizer can minimize any RooAbsReal function with respect to
// its parameters. Usual choices for minimization are RooNLLVar
// and RooChi2Var
// <p>
// RooMinimizer has methods corresponding to MINUIT functions like
// hesse(), migrad(), minos() etc. In each of these function calls
// the state of the MINUIT engine is synchronized with the state
// of the RooFit variables: any change in variables, change
// in the constant status etc is forwarded to MINUIT prior to
// execution of the MINUIT call. Afterwards the RooFit objects
// are resynchronized with the output state of MINUIT: changes
// parameter values, errors are propagated.
// <p>
// Various methods are available to control verbosity, profiling,
// automatic PDF optimization.
// END_HTML
#ifndef __ROOFIT_NOROOMINIMIZER
#include "RooFit.h"
#include "Riostream.h"
#include "TClass.h"
#include <fstream>
#include <iomanip>
#include "TH1.h"
#include "TH2.h"
#include "TMarker.h"
#include "TGraph.h"
#include "TStopwatch.h"
#include "TDirectory.h"
#include "TMatrixDSym.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooAbsReal.h"
#include "RooAbsRealLValue.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooSentinel.h"
#include "RooMsgService.h"
#include "RooPlot.h"
#include "RooMinimizer.h"
#include "RooFitResult.h"
#include "Math/Minimizer.h"
#if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3)
char* operator+( streampos&, char* );
#endif
using namespace std;
ClassImp(RooMinimizer)
;
ROOT::Fit::Fitter *RooMinimizer::_theFitter = 0 ;
void RooMinimizer::cleanup()
{
if (_theFitter) {
delete _theFitter ;
_theFitter =0 ;
}
}
RooMinimizer::RooMinimizer(RooAbsReal& function)
{
RooSentinel::activate() ;
_extV = 0 ;
_func = &function ;
_optConst = kFALSE ;
_verbose = kFALSE ;
_profile = kFALSE ;
_profileStart = kFALSE ;
_printLevel = 1 ;
_minimizerType = "Minuit";
if (_theFitter) delete _theFitter ;
_theFitter = new ROOT::Fit::Fitter;
_fcn = new RooMinimizerFcn(_func,this,_verbose);
_theFitter->Config().SetMinimizer(_minimizerType.c_str());
setEps(1.0);
_theFitter->Config().MinimizerOptions().SetMaxIterations(500*_fcn->NDim());
_theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500*_fcn->NDim());
setPrintLevel(-1) ;
setErrorLevel(_func->defaultErrorLevel()) ;
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
if (RooMsgService::instance().silentMode()) {
setPrintLevel(-1) ;
} else {
setPrintLevel(1) ;
}
}
RooMinimizer::~RooMinimizer()
{
if (_extV) {
delete _extV ;
}
if (_fcn) {
delete _fcn;
}
}
void RooMinimizer::setStrategy(Int_t istrat)
{
_theFitter->Config().MinimizerOptions().SetStrategy(istrat);
}
void RooMinimizer::setMaxIterations(Int_t n)
{
_theFitter->Config().MinimizerOptions().SetMaxIterations(n);
}
void RooMinimizer::setMaxFunctionCalls(Int_t n)
{
_theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(n);
}
void RooMinimizer::setErrorLevel(Double_t level)
{
_theFitter->Config().MinimizerOptions().SetErrorDef(level);
}
void RooMinimizer::setEps(Double_t eps)
{
_theFitter->Config().MinimizerOptions().SetTolerance(eps);
}
void RooMinimizer::setOffsetting(Bool_t flag)
{
_func->enableOffsetting(flag) ;
}
void RooMinimizer::setMinimizerType(const char* type)
{
_minimizerType = type;
}
ROOT::Fit::Fitter* RooMinimizer::fitter()
{
return _theFitter ;
}
const ROOT::Fit::Fitter* RooMinimizer::fitter() const
{
return _theFitter ;
}
RooFitResult* RooMinimizer::fit(const char* options)
{
TString opts(options) ;
opts.ToLower() ;
if (opts.Contains("v")) setVerbose(1) ;
if (opts.Contains("t")) setProfile(1) ;
if (opts.Contains("l")) setLogFile(Form("%s.log",_func->GetName())) ;
if (opts.Contains("c")) optimizeConst(1) ;
if (opts.Contains("0")) setStrategy(0) ;
migrad() ;
if (opts.Contains("0")) setStrategy(1) ;
if (opts.Contains("h")||!opts.Contains("m")) hesse() ;
if (!opts.Contains("m")) minos() ;
return (opts.Contains("r")) ? save() : 0 ;
}
Int_t RooMinimizer::minimize(const char* type, const char* alg)
{
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
_theFitter->Config().SetMinimizer(type,alg);
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
bool ret = _theFitter->FitFCN(*_fcn);
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("MINIMIZE",_status) ;
return _status ;
}
Int_t RooMinimizer::migrad()
{
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad");
bool ret = _theFitter->FitFCN(*_fcn);
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("MIGRAD",_status) ;
return _status ;
}
Int_t RooMinimizer::hesse()
{
if (_theFitter->GetMinimizer()==0) {
coutW(Minimization) << "RooMinimizer::hesse: Error, run Migrad before Hesse!"
<< endl ;
_status = -1;
}
else {
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str());
bool ret = _theFitter->CalculateHessErrors();
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("HESSE",_status) ;
}
return _status ;
}
Int_t RooMinimizer::minos()
{
if (_theFitter->GetMinimizer()==0) {
coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!"
<< endl ;
_status = -1;
}
else {
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str());
bool ret = _theFitter->CalculateMinosErrors();
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("MINOS",_status) ;
}
return _status ;
}
Int_t RooMinimizer::minos(const RooArgSet& minosParamList)
{
if (_theFitter->GetMinimizer()==0) {
coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!"
<< endl ;
_status = -1;
}
else if (minosParamList.getSize()>0) {
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
TIterator* aIter = minosParamList.createIterator() ;
RooAbsArg* arg ;
std::vector<unsigned int> paramInd;
while((arg=(RooAbsArg*)aIter->Next())) {
RooAbsArg* par = _fcn->GetFloatParamList()->find(arg->GetName());
if (par && !par->isConstant()) {
Int_t index = _fcn->GetFloatParamList()->index(par);
paramInd.push_back(index);
}
}
delete aIter ;
if (paramInd.size()) {
_theFitter->Config().SetMinosErrors(paramInd);
_theFitter->Config().SetMinimizer(_minimizerType.c_str());
bool ret = _theFitter->CalculateMinosErrors();
_status = ((ret) ? _theFitter->Result().Status() : -1);
_theFitter->Config().SetMinosErrors(kFALSE);
}
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("MINOS",_status) ;
}
return _status ;
}
Int_t RooMinimizer::seek()
{
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek");
bool ret = _theFitter->FitFCN(*_fcn);
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("SEEK",_status) ;
return _status ;
}
Int_t RooMinimizer::simplex()
{
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex");
bool ret = _theFitter->FitFCN(*_fcn);
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("SEEK",_status) ;
return _status ;
}
Int_t RooMinimizer::improve()
{
_fcn->Synchronize(_theFitter->Config().ParamsSettings(),
_optConst,_verbose) ;
profileStart() ;
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
RooAbsReal::clearEvalErrorLog() ;
_theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved");
bool ret = _theFitter->FitFCN(*_fcn);
_status = ((ret) ? _theFitter->Result().Status() : -1);
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
profileStop() ;
_fcn->BackProp(_theFitter->Result());
saveStatus("IMPROVE",_status) ;
return _status ;
}
Int_t RooMinimizer::setPrintLevel(Int_t newLevel)
{
Int_t ret = _printLevel ;
_theFitter->Config().MinimizerOptions().SetPrintLevel(newLevel+1);
_printLevel = newLevel+1 ;
return ret ;
}
void RooMinimizer::optimizeConst(Int_t flag)
{
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
if (_optConst && !flag){
if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: deactivating const optimization" << endl ;
_func->constOptimizeTestStatistic(RooAbsArg::DeActivate) ;
_optConst = flag ;
} else if (!_optConst && flag) {
if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: activating const optimization" << endl ;
_func->constOptimizeTestStatistic(RooAbsArg::Activate,flag>1) ;
_optConst = flag ;
} else if (_optConst && flag) {
if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: const optimization already active" << endl ;
} else {
if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: const optimization wasn't active" << endl ;
}
RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
}
RooFitResult* RooMinimizer::save(const char* userName, const char* userTitle)
{
if (_theFitter->GetMinimizer()==0) {
coutW(Minimization) << "RooMinimizer::save: Error, run minimization before!"
<< endl ;
return 0;
}
TString name,title ;
name = userName ? userName : Form("%s", _func->GetName()) ;
title = userTitle ? userTitle : Form("%s", _func->GetTitle()) ;
RooFitResult* fitRes = new RooFitResult(name,title) ;
Int_t i ;
RooArgList saveConstList(*(_fcn->GetConstParamList())) ;
RooArgList saveFloatInitList(*(_fcn->GetInitFloatParamList())) ;
RooArgList saveFloatFinalList(*(_fcn->GetFloatParamList())) ;
for (i=0 ; i<_fcn->GetFloatParamList()->getSize() ; i++) {
RooAbsArg* par = _fcn->GetFloatParamList()->at(i) ;
if (par->isConstant()) {
saveFloatInitList.remove(*saveFloatInitList.find(par->GetName()),kTRUE) ;
saveFloatFinalList.remove(*par) ;
saveConstList.add(*par) ;
}
}
saveConstList.sort() ;
fitRes->setConstParList(saveConstList) ;
fitRes->setInitParList(saveFloatInitList) ;
fitRes->setStatus(_status) ;
fitRes->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ;
fitRes->setMinNLL(_theFitter->Result().MinFcnValue()) ;
fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()) ;
fitRes->setEDM(_theFitter->Result().Edm()) ;
fitRes->setFinalParList(saveFloatFinalList) ;
if (!_extV) {
std::vector<double> globalCC;
TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ;
TMatrixDSym covs(_theFitter->Result().Parameters().size()) ;
for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) {
globalCC.push_back(_theFitter->Result().GlobalCC(ic));
for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) {
corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii);
covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii);
}
}
fitRes->fillCorrMatrix(globalCC,corrs,covs) ;
} else {
fitRes->setCovarianceMatrix(*_extV) ;
}
fitRes->setStatusHistory(_statusHistory) ;
return fitRes ;
}
RooPlot* RooMinimizer::contour(RooRealVar& var1, RooRealVar& var2,
Double_t n1, Double_t n2, Double_t n3,
Double_t n4, Double_t n5, Double_t n6)
{
RooArgList* params = _fcn->GetFloatParamList() ;
RooArgList* paramSave = (RooArgList*) params->snapshot() ;
Int_t index1= _fcn->GetFloatParamList()->index(&var1);
if(index1 < 0) {
coutE(Minimization) << "RooMinimizer::contour(" << GetName()
<< ") ERROR: " << var1.GetName()
<< " is not a floating parameter of "
<< _func->GetName() << endl ;
return 0;
}
Int_t index2= _fcn->GetFloatParamList()->index(&var2);
if(index2 < 0) {
coutE(Minimization) << "RooMinimizer::contour(" << GetName()
<< ") ERROR: " << var2.GetName()
<< " is not a floating parameter of PDF "
<< _func->GetName() << endl ;
return 0;
}
RooPlot* frame = new RooPlot(var1,var2) ;
TMarker *point= new TMarker(var1.getVal(), var2.getVal(), 8);
frame->addObject(point) ;
Double_t errdef= _theFitter->GetMinimizer()->ErrorDef();
Double_t n[6] ;
n[0] = n1 ; n[1] = n2 ; n[2] = n3 ; n[3] = n4 ; n[4] = n5 ; n[5] = n6 ;
unsigned int npoints(50);
for (Int_t ic = 0 ; ic<6 ; ic++) {
if(n[ic] > 0) {
_theFitter->GetMinimizer()->SetErrorDef(n[ic]*n[ic]*errdef);
Double_t *xcoor = new Double_t[npoints+1];
Double_t *ycoor = new Double_t[npoints+1];
bool ret = _theFitter->GetMinimizer()->Contour(index1,index2,npoints,xcoor,ycoor);
if (!ret) {
coutE(Minimization) << "RooMinimizer::contour("
<< GetName()
<< ") ERROR: MINUIT did not return a contour graph for n="
<< n[ic] << endl ;
} else {
xcoor[npoints] = xcoor[0];
ycoor[npoints] = ycoor[0];
TGraph* graph = new TGraph(npoints+1,xcoor,ycoor);
graph->SetName(Form("contour_%s_n%f",_func->GetName(),n[ic])) ;
graph->SetLineStyle(ic+1) ;
graph->SetLineWidth(2) ;
graph->SetLineColor(kBlue) ;
frame->addObject(graph,"L") ;
}
delete [] xcoor;
delete [] ycoor;
}
}
_theFitter->Config().MinimizerOptions().SetErrorDef(errdef);
*params = *paramSave ;
delete paramSave ;
return frame ;
}
void RooMinimizer::profileStart()
{
if (_profile) {
_timer.Start() ;
_cumulTimer.Start(_profileStart?kFALSE:kTRUE) ;
_profileStart = kTRUE ;
}
}
void RooMinimizer::profileStop()
{
if (_profile) {
_timer.Stop() ;
_cumulTimer.Stop() ;
coutI(Minimization) << "Command timer: " ; _timer.Print() ;
coutI(Minimization) << "Session timer: " ; _cumulTimer.Print() ;
}
}
void RooMinimizer::applyCovarianceMatrix(TMatrixDSym& V)
{
_extV = (TMatrixDSym*) V.Clone() ;
_fcn->ApplyCovarianceMatrix(*_extV);
}
RooFitResult* RooMinimizer::lastMinuitFit(const RooArgList& varList)
{
if (_theFitter==0 || _theFitter->GetMinimizer()==0) {
oocoutE((TObject*)0,InputArguments) << "RooMinimizer::save: Error, run minimization before!"
<< endl ;
return 0;
}
if (varList.getSize()>0 && varList.getSize()!=Int_t(_theFitter->Result().NTotalParameters())) {
oocoutE((TObject*)0,InputArguments)
<< "RooMinimizer::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl
<< " or match the number of variables of the last fit ("
<< _theFitter->Result().NTotalParameters() << ")" << endl ;
return 0 ;
}
TIterator* iter = varList.createIterator() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)iter->Next())) {
if (!dynamic_cast<RooRealVar*>(arg)) {
oocoutE((TObject*)0,InputArguments) << "RooMinimizer::lastMinuitFit: ERROR: variable '"
<< arg->GetName() << "' is not of type RooRealVar" << endl ;
return 0 ;
}
}
delete iter ;
RooFitResult* res = new RooFitResult("lastMinuitFit","Last MINUIT fit") ;
RooArgList constPars("constPars") ;
RooArgList floatPars("floatPars") ;
UInt_t i ;
for (i = 0; i < _theFitter->Result().NTotalParameters(); ++i) {
TString varName(_theFitter->Result().GetParameterName(i));
Bool_t isConst(_theFitter->Result().IsParameterFixed(i)) ;
Double_t xlo = _theFitter->Config().ParSettings(i).LowerLimit();
Double_t xhi = _theFitter->Config().ParSettings(i).UpperLimit();
Double_t xerr = _theFitter->Result().Error(i);
Double_t xval = _theFitter->Result().Value(i);
RooRealVar* var ;
if (varList.getSize()==0) {
if ((xlo<xhi) && !isConst) {
var = new RooRealVar(varName,varName,xval,xlo,xhi) ;
} else {
var = new RooRealVar(varName,varName,xval) ;
}
var->setConstant(isConst) ;
} else {
var = (RooRealVar*) varList.at(i)->Clone() ;
var->setConstant(isConst) ;
var->setVal(xval) ;
if (xlo<xhi) {
var->setRange(xlo,xhi) ;
}
if (varName.CompareTo(var->GetName())) {
oocoutI((TObject*)0,Eval) << "RooMinimizer::lastMinuitFit: fit parameter '" << varName
<< "' stored in variable '" << var->GetName() << "'" << endl ;
}
}
if (isConst) {
constPars.addOwned(*var) ;
} else {
var->setError(xerr) ;
floatPars.addOwned(*var) ;
}
}
res->setConstParList(constPars) ;
res->setInitParList(floatPars) ;
res->setFinalParList(floatPars) ;
res->setMinNLL(_theFitter->Result().MinFcnValue()) ;
res->setEDM(_theFitter->Result().Edm()) ;
res->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ;
res->setStatus(_theFitter->Result().Status()) ;
std::vector<double> globalCC;
TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ;
TMatrixDSym covs(_theFitter->Result().Parameters().size()) ;
for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) {
globalCC.push_back(_theFitter->Result().GlobalCC(ic));
for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) {
corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii);
covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii);
}
}
res->fillCorrMatrix(globalCC,corrs,covs) ;
return res;
}
#endif